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' ABSTRACT 


Supersonic flow past oscillating flat plate cascades 
with supersonic leading-edge locus is analysed using a lin- 
earized method of characteristics valid for arbitrary fre- 
quencies and an elementary analytical theory valid only for 
low frequencies of oscillation. These two methods are ex- 
tensions of previous work by Teipel and Sauer for the single 
airfoil in an unbounded supersonic flow to the case of air- 
foils oscillating in cascade. Included is the determination 
of pressure distribution and both a two-degree-of-freedom 
(bending and torsion) Elmo ces analysis and a single-degree- 
of-freedom (torsion) flutter analysis. Numerically deter- 
mined flutter boundaries —_ presented for various primary 
parameters onic as, Mach number, solidity, stagger angle, 
density ratio, structural damping coefficient, and elastic 
axis position. Also, results are presented for the related 
problem of supersonic wind tunnel interference (including 
the effect of tunnel porosity) and airfoil-airfoil interfer- 


ence. 
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L. INTRODUCTION 


In order to design an effective supersonic compressor, an 
analysis of the complex unsteady flow phenomena in this speed 
regime is necessary. Since the infinite cascade has long been 
used by the engineer to model two-dimensional compressor 
flows, the study of oscillating supersonic cascades is import- 
ant in predicting the flutter characteristics or dynamic res- 
ponse of these compressors. Moreover, the basic understanding 
of the flow gained by such a study can lead to better design 
criteria. 

Two basic cascade configurations can be distinguished at 
the outset; i.e., cascades with either subsonic or supersonic 
leading-edge locus. Although the former is the case of pri- 
Mary interest in current research, in this paper the case of 
Supersonic leading-edge locus is treated in ite 6 form a 
basis for extension to the more complicated case of subsonic 
leading-edge locus. This method of approach also has the 
advantage that the problem of oscillatory wind tunnel inter- 
ference, already analysed in previous work, is contained in 
the theory as a special case. 

The interference problem of linearized supersonic flow 
past airfoils oscillating between solid wind tunnel walls 
was considered by Miles (1956) who derived a solution using 
Laplace transform techniques. Drake (1956) treated this case 
for wind tunnels with free jet boundaries and later Drake 


(1957) gave a solution to this interference problem for 





porous-wWall wind tunnels also using Laplace transform methods. 
In an extension of Miles’ work, Lane (1957) presented a solu- 
tion for supersonic flow past oscillating cascades with super- 
sonic leading-edge locus again using Laplace transform tech- 
niques. Further work in this area was done by Hamamoto (1962). 
Using Teipel's (1962) method of characteristics approach, 
Platzer and Pierce (1970) made an analysis of oscillatory 
Supersonic wind tunnel interference and with the help of the. 
high speed computer were able to predict the pressure distri- 
bution along an airfoil oscillating with arbitrary frequency 
in solid, free jet, or porous-wall wind tunnels. Platzer 
(1971) will present an elementary analysis of porous-wall 
wind tunnel interference effects which generalizes Sauer's 
(1950) solution for a single airfoil oscillating at low fre- 
quency in an unbounded linearized supersonic flow. Platzer 
and Chalkley (1972) further extended this solution to form an 
elementary theory that together with a method of characteris- 
‘tics procedure (based on Platzer and Pierce, 1970) is used to 
analyse the supersonic flow past oscillating flat plate cas- 
cades having supersonic axial velocities but otherwise arbi- 
trary stagger angle. 

In this thesis the detailed description of both the ele- 
mentary theory and the method of characteristics that was not 
possible in Platzer and Chalkley (1972) is given. With the 
aid of the high speed computer, the method of characteristics 
procedure is used to predict the pressure distributions and 


aerodynamic forces and pitching moments on a typical cascade 





blade. Both a two-degree-of-freedom (torsion and bending) 
flutter analysis and a one-degree-of-freedom (torsion) flut- 
ter analysis are presented, with numerical results given for 
the case of the torsional flutter. The elementary theory 
(valid only for low frequencies of oscillation) is used to 
provide a convenient check of the characteristics approach 
as well as to easily predict a number of dynamic instability 
boundaries for a slowly oscillating airfoil subjected to 
interference from supersonic wind tunnel walls in one case 
and from a larger airfoil in close proximity in yet another 


case. 
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IL. PROBLEM FORMULATION 


The cascade is considered to be one of two-dimensional 
flat plates with each blade performing equal low amplitude, 
Simple harmonic oscillations of the same mode with the same 
interblade phase angle. As shown in Figure 1, it is aligned 
in the x-y coordinate system such that one blade lies along 


the x-axis with its leading-edge at the origin. 
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The cascade is further considered to have supersonic 
leading-edge locus and, except for this restriction, arbi- 
trary stagger, such that, 

tan BS cot a (1) 
With this assumption, the flow between any two adjacent 
blades can be used to describe the flow through the entire 


cascade. Since disturbances cannot travel upstream of the 


eh 








Mach lines from the leading-edge of a blade, interference 
from other than adjacent paeaee is not possible in the super- 
sonic flow between them. Hence, only the flow between the 
blade along the x-axis and the one above it is considered. 
This flow is assumed to be the non-viscous flow of a perfect 


gas governed by the continuity equation, 


Dp Cachan 


the Euler equations, 


ee ee 
Dt oS 5 OR 0 (3a) 
UNE th i 
Teo os 0 (3b) 


and the energy equation, 


Ds 
Dt 


=a (4) 
Because of the assumption of small amplitude oscillations 
of the cascade blades, all flow quantities are considered to 
be small perturbations linearly superimposed on freestream 
quantities. The velocities are written as 
u=u tu (5a) 
c=e te (5b) 
while the pressure and density perturbations are written as 


ye Pp -— Pz (6a) 


MPs 9 - Py, (6b) 


With the assumption of isentropic flow, the local velo- 


city of sound can be defined as, 


we (7) 


az 








har ther, 


—F_ «= constant (8) 
oY 
Taking the total differential, 
~ i dp - yp 1. do = 0 C9) 
pi ie 
or 
Py Pp 
> (10) 
hence, 
Coma y a (11) 


The flow boundary condition along the blade surfaces [as 
Shown in Ch. 5, Bisplinghoff, Ashley, Halfman (1955) fora 
simple airfoil] is as follows: If the surface over which the 
fluid flows has the equation, 

Fase y, t) = 0 (12) 
particles in contact with the surface must have the same nor- 
mal velocity as the surface. Stated differently, the rate of 
change of F is zero when the motion of a particular fluid 
element is followed along the surface, or 


DF _ 
DL = Q (13) 


If an arbitrary, thin, two-dimensional cascade blade is 
considered, the equation of the upper surface can be written, 


F Gx yt) ay Se (x,t) = 0 (14) 


maere y = 0 is the mean camber line and Ya is the distance 
from the mean camber line to the upper surface. Likewise, 


the equation of the lower surface can be written, 


eS 
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At y = ve 
en ie 
Dt a 
and at y = Yy 
ma. ii 
Dt ot 
since 
ig 
dy dy 
Thus, the normal flow vel 
ot ox 
dt ox 


By applying 


freestream, the normal ve 
dt co gx 
dt co gx 

Since a thin blade is bei 


thus the equations become, 


dt o ox 
oes 
dt o ot 


This normal yelocity 


Series expansion of the normal flow velocity at y 


= y¥ - yt 0 Gis) 
5 Oe 
u Pa + y= QO (16) 
” 9Yy 
eee = 0 (17) 

(18) 

ocity can be written, 
at y = y,6*st) (19a) 
at oy = y, (x,t) (73ib») 


the assumption of linear perturbations on the 


locity equations become 


ay. 
Pg at y = y, (x,t) (20a) 
oy, 
ok at y = y, (x,t) (20b) 


ng considered vie and y, are small, 


cancelling the higher order terms, 
at y = y,, 6%st) (21a) 
at y = y, («5 t) (21b) 


can be further written as a Taylor 


QO. Thus, 
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peer 
GLO ey + 9 onalee Os, t) 


v(x,y >t) = dy 
ied sauessaiis 
to7 i Ht ee C223) 
: Jy? 
and | 
- dv(x,0 ,t) 
= ee ee 
v(x,y,5t) v(x,0 ,t) +t y, By 
a ah A © FD, 
t oT dy? ae eoeeeee (22b) 


Again with the assumption of a thin blade and small linear 


perturbations on the freestream: y_, Yy> and v are small and 


u 


higher order terms may again be cancelled from the equations 


leaving: 


v(x,y 2) - v(x,0',t) (23a) 
v(x,y,,t) = v(x,0 ,t) (23b) 
Thus, 
dy 
v(x,y,t) = ee + uy = at y = aad (24a) 
oY, oy, 7 
v(x,y,t) = we ea at y = 0 (24b) 


With the assumption of simple harmonic motion of the blade, 


y 0) aoe (25a) 


y 6% t) 


LWwt 


y, (x,t) y, (x) e C255) 


and the normal velocity equations can then be written, 


dy 

v(x,y,t) = Liwy () fae = ane at y = on (26a) 
OY, iwt — 

v(x,y,t) = [iwy, (x) ieee ree aa at y = 0 (26b) 


Since the flow is identical between any two adjacent 


blades, 
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See er L a 
v(x,y,t) = fiwy, (x) he 7x Je ate y= .d C2) 
In plunge, the lower blade oscillation about y = 0 is, 
y= - hy et¥F (28) 
hence, 
v(x,y,t) = [-iwh ] e*¥F at y= 0° (29) 
while the upper blade oscillation about y = dis, 
ae ho ot (wtts) (30) 


where § is the interblade phase angle. Hence, 


v(x,y,t) = [w sin 4 ho - iw cos 64 hol a7 UE (31) 
at y = d- 
In pitch, the lower blade oscillation about y = 0 is, 
zh iwt 
yr- Q | (x- x) e (32) 
hence, 
te = [wee = tue Gees (33) 
aed co "O° O O 


at y = Oe: while the upper blade oscillation about y = d is, 


y= -8 (x - B-x,) e \Wtto) (34) 
hence, 
v(x,y,t) = ae {um cos 6 -wW sin 6 Cee 
- 16) cue sin 6 + Ww cos 6 (x-B-x 9) }] Bee C33) 
at y = d- 


The pressure distribution on the blades can be written in 
terms of the local sonic velocity as follows: In terms of 


perturbation quantities, Eq. (11) can be written, 
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1 P, + Ap 


Po (1 + ve, 


2 = 
Santezcc. = Y 


oO 


(36) 


cancelling the higher order c* term. Since Be is less than 


co 


one, Eq. (36) can be expanded in the geometric series, 


‘ 2 
4 - a Oe - (37) 
(1+ =P) ea Pe Poo 
Peo 


neglecting the higher order terms. Substituting this in 


Eq. (36) gives, 








1 © Ap 
2 = = lS A 38 
Since Ap and Ap are small, 
4p, . 1 
Ap” dp a2 7 — 
dp 


Substituting this in Eq. (38) gives after some algebraic 
manipulation, 


1 1 





/ _ 
yo1 cc. S me: ) = >. Ap (40) 
Cayce. 
Since oy ae ye is small, st - may be expanded in 
So 
(y-1)c, 
the geometric series, 

1 2c 2c 2 7a (03 
—_—— -hOCTaS cL hl r CFO a fae ae 6 ee 6) ae ine Ce 
A, .2e (y-ie, * 'Qy=tyer! ae G aT) 

Sle 
(41) 


neglecting higher order terms. Hence, Eq. (40) can be writ- 


ten, neglecting higher order terms, 


Pacer 7 2, 'c, (c - cz) (42) 
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Finally, any non-dimensional quantities so stated in the 
paper are made such with reference to the blade chord length 


and the freestream velocity unless otherwise defined. 
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Ifl. METHOD OF CHARACTERISTICS 


In the method of characteristics it is desirable to as- 
certain the coordinate system across which all possible dis- 
continuities in flow properties may occur. Along this 
coordinate system the equations of moeion of the flow field 
can then be treated as ordinary differential equations that 
are solvable by classical or numerical techniques (e.g., 
finite differences). To obtain the equations of this coor- 
dinate system (the characteristic directions) in the (x,y) 
plane, the equations of motion are written in terms of the 


arbitrary intersecting coordinates, 


Sg REISS), (43a) 


and 


n= 1 55) (43b) 
If the first derivatives of the dependent variables, u, Mas 
and p (or c) with respect to & are made indeterminate across 
lines of n = constant, and the first derivatives of the de- 
pendent variables with respect to n are made indeterminate 
across lines of & = constant, then any possible discontinui- 
ties in the first derivatives will occur across these lines. 
These lines are then the characteristics and their equations 
are obtained in the evaluation of the indeterminacies. 


Consider, first, two-dimensional, steady flow. The gov- 


erning equations of motion are the continuity equation, 


| ch — =O (44) 


uo 





the Euler equations, 


— ou du, 1 ap _ 
u a5 + v ay + 5 OES L 0 (45a) 
- ov, av, 1 dp _ 
u 35 + ¥ ay + ae 0 (45b) 
and the energy equation, 
- ds 0S _ 
ua ay 55) 0 (46) 
Along a streamline this equation may be written, 
u SP + y SP cya 22 + y (47) 


ox ay oy 


Substituting in Eq. (44), the continuity equation becomes, 


oN cee tee sy, ORG (48) 


as oy Ae dX dy 


OD 
== 
ates 
~lc! 


In terms of the arbitrary coordinates, Eqs. (43), the 


continuity equation becomes, 


= 1 = 
pe Ur a8 12. Sa 15 = [ug fe viene Mls 


1 
peewee, Um PY. Ve ree 
Cc 


ae Sel [u n, tv ed Pp (49) 


2 ia 


In like manner, the Euler equations become, 


pug. a WS el EP ry = -p[un, a vn] ae Peay (50a) 


S is) 


pfug, + ele y = Sele aa Ti SLi (50b) 


ge Ele ye 


Eq. (49) and Eqs. (50) form a system of three simultaneous 


.§ 


rule gives a ratio of two determinants. The denominator is 


equations in Me? v_, and Pe Solving for Pe by Cramer's 


the determinant of the matrix of coefficients of Ur, Ves and 


Pr formed by Eq. (49) and Eqs. (50), while the numerator is 
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this determinant with the right hand side of the equations 
substituted for the coefficients of Pe 
Since across lines of yn = constant, Pe is indeterminate, 
both the numerator and the denominator must equal zero. Set- 
ting the denominator equal to zero and expanding the deter- 
minant gives, 

Pate. + vey] [(u*- c7) ES + 2uv bes +(v*- e*)Er] = 0 

(2) 

The solutions to this equation give the equations of all 


three characteristics in the physical (x,y) plane. These are, 


Cre 
7 ae tan z | (52a) 
a = tan (z +t a) (52b) 


where zg is the angle the streamline makes with the x-axis and 
@ is the Mach angle. All three characteristic directions are 
obtained from Eq. (51) because &€ and Nn are arbitrary and in- 
terchangeable: If Eq. (49) and Eqs. (50) were written such 
that derivatives with respect to n are the unknowns, making 
Py indeterminate across lines of — = constant would result in 
msarutcions identical to Eqs. (52), a priori. Eq. (52a) is 
merely the equation of a streamline while Eqs. (52b) are the 
equations of the right and left-running Mach lines. Tradi- 
tionally, Eqs. (52b) are the equations of € = constant and 
iie= COnStant. 

The compatibility relation which relates changes in Ae aa 


and p along the (&,n) characteristics can be obtained by 


Za 





setting the numerator of Py equal to zero. Expanding the 
determinant and solving the resultant equation yields after 


some algebraic manipulation, 


ou ht - l dp 
Se — Seis + —_ — 
Vv 5 u cot @ 5 0 (53) 


This derivation is shown in detail in Chapter 8 of Oswatitsch 
(1956). These relations must be satisfied along the (&,n) 
Gnarwacteristics. 

The physical significance of these characteristics is 
that they act as boundaries across which any discontinuities 
will occur. In this flow field where streamlines and Mach 
lines are the characteristics, a discontinuity is readily 
apparent across the Mach line emanating from the leading-edge 
of a slender airfoil or across a slip plane (streamline) 
emanating from a Mach stem. Characteristics can further be 
described as information carriers in that along them all dis- 
turbances propagate and all flow quantities can be determined 
given their values in the domain of dependence. It should 
also be noted that the characteristic directions are com- 
pletely independent of the type of coordinate system used to 
describe the flow. They are determined solely by the physical 
nature of the flow field as expressed in the equations of 
motion. 

The problem of two-dimensional, unsteady flow over a flat 
plate as shown by Teipel (1962) can be considered in the same 
Manner as before. The governing equations of motion [Eq. (2), 


Eq. (3), and Eq. (4)] are written in an arbitrary coordinate 
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System across which variation of flow properties are indeter- 
minate. The evaluation of these indeterminacies will again 
give the characteristic directions of these coordinates. 

The: arbitrary coordinates are given in Eqs. (48). Eq.(2), 
Eq. (3), and Eq. (4) can be rewritten in a form more suitable 
for transformation by using the local sonic velocity as an 
unknown. Assuming that air is a perfect gas, Eq. (11) can 
be written, | 

c* = ¥ & it (54) 
where @® is the universal gas constant and is the molecular 
weight. Taking the total differential and dividing by c* 
gives, 
dc R aT 


2 = Va 25 (55) 


Q 


Substituting Eq. (54) into this gives, 


C dT 
T (56) 


In like manner, the total differential of Eq. (11) is, 


2ede = y BE (57) 


mad £rom Eq. (11) again, 


(58) 


fo 
Bt bs 
O1 
1! 
fe 
i 
hay 


From the First Law of Thermodynamics and the definition of 
entropy, 


Tds = du -FP do (59) 
(elim 
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Since du = * dt,. this may be written as, (60) 


4 Be ko 
Tds c dT RT ; (61) 


where R is the gas constant Cc ane and or and care the 
specific heats of the gas at constant pressure and volume, 
Bespectiyvely. Dividing Eq. (61) by Te, substituting from 


Eq. (56) and Eq. (59), and then taking the substantial deri- 





vative of S with respect to time gives, 


Vv 


QO} 


D 


1 m 
ra = = Dt 


Vv 


LIX 
rin 
Ont 


2 Do 
Se aye = (62) 


since entropy is held constant. This can be written substi- 


mmtane from Eq. (11), 


D 2 ~ De 
Be. 2 oe BE (63) 


Introducing the small perturbation assumption and substitut- 
ing from Eq. (7) and Eq. (63), Eq. (2) can be written, can- 


celling higher order terms, 


om 


3 Z dc 4) (64) 


ot. 





Cc 
x | 

ane 

Q 
male 
“iC 

te 

n 
ele 
<1< 


2 
y~-1 


Q 
ct 
1 
_— 
8 


This is the form of the continuity equation that will be used 
in the coordinate transformation. Again dividing Eq. (61) by 
Te, substituting from Eq. (56) and Eq. (58), and then taking 
the partial derivative with respect to x (holding entropy 


constant) gives, 


M1 


ae. 2 32 


eee 
Lo 


ox y-1 


| 


DlK 


(65) 


Qz 


x 
Substituting in the first Euler equation, Eq. (3) becomes, 


cancelling higher order terms, 
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du du, 2 Oe. 
aE se thes = 1 ET Caria 0 (66) 


In like manner, the second Euler equation, Eq. (4), becomes, 





c ww = 0 (67) 


The continuity equation and the Euler equations thus form a 
system of three simultaneous equations in u, VV, and Ce Wen 
assumption of simple harmonic motion, Teipel (1962) intro- 


duced his amplitude functions, 











iwt Ee Mo 
U(x,ye OF = ——— (68) 
He po ————_— = (69) 
YM*-] Yeo 
CC eae = 2 lees ce (70) 
¥=1 <4,25 ce 


where U, V, and C are complex. 


Substituting in Eq. (64), Eq. (66), and Eq. (67) gives, 





ou y INGE Ce ee ro Sa Cia) 
OX dy x 
OU pee. ane (72) 
Ox OX 
ov, 1 CA (73) 


Ox /M2-1 oy 
where k is the reduced frequency. 
Transforming this system into the new coordinate systen, 
E(x,y) and n(x,y), and considering now a system of three 


simultaneous equations in unknowns, Dee Ufa and ee the equa- 


tions may be written, (74) 
2 2 a: 
= -U_- t a -M“C - ikM’C 
a U, + cot a Doel faa ME Se n co Vany nex 
— = =— - i U Vis 
ex U; + Ee Cy cra cin, ik (753 
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Y_ttan a C= - _— . @ — key, 
bx 2 ee Sy ¢ 42 ny tam @ C Ny iky, (7B9. 


As before Urs Ves and oe are made indeterminat® arross:; 
dines of n = constant. This forces any discontimitiess im 
U, V, and C to lie across Lines of n = constant. Ome condi-- 
tion of indeterminacy is that the determinant of the: matrix 
of coefficients of Urs Ve» and Cr in the above equatcionss: be: 
zero. If this determinant is expanded and then set eyurad to: 
zero, the resulting equation is, 


Sec ie, ete = 0 (97) 


x 
As before the solutions to this equation give the thre#: char>- 


-acteristic directions in the physical (x,y) plane: 


d = ~ = 
ae 0 (7Ba) 
d =_— ’ = 
Ser + tan a (7Bb:). 


These coincide with the results af the steady flow. analy-- 
sis when it is remembered that the small perturbation assump-- 
tion forces the streamlines to be parallel with the x-axis 
(c = 0). 

In order to find the compatibility relations of the first 


derivatives of U, V, and C along the characteristics, the 


equations of motion are written along lines of € = constant, 
N = constant, and along streamlines, that is, for 
q 
str(x,y) = constant, == 0 (79a) 
e (x,y) = constant, ages ee (79b,) 


Bib 








ncx,y) = constant, $2 = sees (79c) 


¢ YM*-1 
If an arbitrary function of x and y in the equations is 


denoted by £ then, ces holding str constant becomes, 





Ox , 
le Oe 
Soe? Ste ox oP) 
Likewise, ved holding & constant becomes, 
= = of 4 i for (81) 
xe ame ak tet ed 
Of , 
and, oa holding n constant becomes, 
of Jf Z. Of 
-, See Se 82 
Ox 1 ox J/M2-1 oy ca 
Combining these equations and rewriting, 
of 2 ae 
ox (Ox str Se) 
OE ale af 
Be peo ee ae 
"7 1, See (83c) 


Substituting these relations in the equations of motion gives, 


1 oU OU 1 OV OV 
21x): + Ga a! * wae = 1) Lay): ~ ee ae 
1 oC oC ; x 
a ee ‘ + Car + ikM*c = 0 (84a) 
1 ,,3U 3U 1,,9¢ 3C ae te 
> [(se + Sree) al + Oe Gaels a Gas) + ikU = 0 (84b) 
OV oy 1 oC aC : _ 
Ga) ai (asda! tf FE ae) is - Care + iky = 0 (84c) 


Subtracting Eq. (84b) from Eq. (84a) gives an equation that 
is first added to Eq. (84c) and then subtracted from Eq. (84c). 


‘The result is two compatibility relations that must be 
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Smevouted along lines of constant — and lines of constant fn. 


These are, a 
i 








(ov. oC 7a 7 

Cris 5 Gee + ik [V + ca (M“C U)] 0 (85a) 
a 

ope = ce a eek Tv y2od (M“C U)] 0 (85b) 


The third compatibility relation is obtained by considering 
Eq. (84b) along a streamline, 


3U rye 
So) eee Ge ey 


I 
-) 


+ ikU (86) 


If irrotationality is assumed, the equation of irrotationali- 
EY, 


Ou OV. 
dy 7 ox ? | (87) 


may be substituted for the second Euler equation in the sys- 
tem. Again with the assumption of simple harmonic motion, it 


may be written in terms of Teipel's amplitude functions, 


le 
! 
2 
j 
_ 
|< 
I 
© 


(88) 


The non-dimensional system of equations can then be written 





as, 
Oe Viel 1 Se oe SS eee ep (89a) 
x y x 
= + oe 4 EE = (89b) 
- ss J/M2 — 1 sv 0 (89c) 


By transforming Eqs. (89) to the arbitrary coordinate systen, 
— and fq, as before, the system can be considered one of simul- 


Caneous equations of U,_, V_, and oo: - By making any of these 


g° 
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meniyativyes indeterminate the results of setting the deter- 
Minant of the matrix of their coefficients equal to zero are 
Eqs. (78), the characteristic directions. 

By writing the equations of motion, as before, along 
lines of € = constant and nN = constant and performing the 
same algebraic manipulations, the two compatibility relations 


can be determined: 








7 ing 
Ge) - G+ ik S W- cc) = 0 (90a) 
5 E M*-1 
gu OV ale 7 
oe a ata =, + ik (CU Cc) = 0 (90b) 
M“=-1 
The third compatibility relation is, as before, 
dU dC : _ 
Som str Cw eer eee SE 


The equations of motion are now reduced to a system con- 
taining only derivatives with respect to x. It is, thus, 
possible to solve these equations using finite Ao kserences 
in the following manner: Separating the real from the imagi- 


nary components, the equations can be written, 


dU, aC, 
oe ie a ome str - NS 8) (91a) 
aU, ocr 
a Str as (Sz SEG os KUR ait (91b) 
dU oV 2 
__& R = M _ ms 
OX Gee ane CUy- Cy) =a SoS) 
M“=-1 
dU oV 2 
L I M 
c-) »- (——) +k CUSEaCe) = 0 (91d) 
oU e) ee. 
R R M 3 = 
ern G ye k a (Uy Cy) = 0 (91le) 


OBS, 





gU, avy M2 
Gen + Te a +k a (U,- C.) =o (91f£) 


To write these equations in finite difference form, a 
computational molecule as shown in Figure 2 will be used. 


P and P known, all the 


With the velocity values at P 12? 21 


ls 


values at Po can be calculated. The distance Ax is arbi- 


Erary. 





Figure 2 


If F is an arbitrary flow quantity, its partial deriva- 


tives can thus be written, 


OF Z ZZ 11 
Ge ace = Ax (92a) 
F - F 
(s= a pee eZ (92b) 
Xx 1 
E ZAx 
F - F 
ES = ai ee (92c) 
ox’ n Ly 
Bia 


The values of the flow quantities in the equations are aver- 
eees, such that, 


(F + F (93a) 


IL3l 22) 
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1 
(F) = 2 (Fi > + F,5) (93b) 
mk 
CF) 3 py + Fo2) (93c) 
With these relations substituted into the equations, the sys- 
tem becomes, after some algebraic manipulation, 
OR “oR eee To R (94a) 
Usor * Saat * ApUooR 7 Siar (94b) 
Usop ~ Yoon ~ ByUaor + Byooat= Kayp 694°) 
Usor ~ Yaar * ByUoop ~ Brooor™ K3gzr 6944) 
Voor + Voor ~ BrUoaz + Bron0r= Ksgp 6948) 
Uoor * Voor * BY oon ~ Brooor Bs5gy ‘6? **? 
where: 
SS eRe alas alibi) ae jal (95a) 
“Wo gp oviase § nae > ora (95b) 
Ryan = Yaar ~ Yuen * Bz ar- Sya7)? 695°) 
Sovage @ stage > ge = By (Uj ope 12R? (95d) 
Roer ~ Yair * Voir * By Your- Caq7) 6958) 
Roex ~ Your * Your ~ Br Uoip™ Soap)? 699F) 
and, 
A. = ak Ax (96a) 
I 2 
ip SONG 
By mei k Ax (96b) 
Men 


As shown in Teipel (1962), 


BL 





solving these equations gives, 





(1-A, 8B, ) [5 (K au +28 


a 


CK 


J 


e.- 34R'56R That L2 asi" BIR 
(1-A,B,)? + (2B,) “ie 

2. — ee WSaviae =361°7*s Kior!7 2b aent tSiar! 
(1-A,B,) + (2B,) (97b) 

Te en (97¢) 

22RD Sone 34K 

Loon = 5 Sam Sayan? SS 

Coot Re 2a ore S2dS) 

CGO «Ire 20 eo EY 


These are the finite difference equations for a point in the 
general flow field. 
At a point on the top of a cascade blade or solid bound- 


ary, the normal velocity V is prescribed by the movement 


Pigg 


of the blade. This velocity was given in the Problem Formu- 


lation. Here the computational molecule is shown in Figure 3. 





Fligaxne> 3 


32 





The applicable equations are, 


On eon «1 22tee OR SY, 
pr ) 20y OR © aan CEP) 
Usorn ~ ByUa07 * BySna7 = K5er~ *34n 6780) 
Voor + ByUgoR ~ BySn2R = Ks6r7 Baar = 6984) 
where Rio and Ree are as before, but, 
oan 7 ae lene 
Sat ai ee) 


as given by the boundary conditions. Solving for these 


equations gives, 


a eta ser 834 or lot 1 eee sone eae eRe 
22R 2 2 
(1 A,BL) + (2B) (100a) 
u 7 (1-A, BL) [Ke 67K yp FB ky op] -2By [Ks 6p -K 3, p7 By Ky 97] 
P21 2 2 ‘ 
efR 15) 2) (100b) 
Poor “12R ~ "20er * 81% 201 (100c) 
Coot Kiar ~ Yoat ~ ArY2er eee? 


Voor and Voor are known from the boundary conditions. These 
are the finite difference equations for the flow velocities 
at the top of a solid boundary or cascade blade. 

Likewise, at the bottom of a solid boundary or cascade 
blade, the normal velocity, Vo9 1s prescribed by the movement 


of the blade, and was given in the Problem Formulation. The 
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computational molecule, different from that at the top of a 


blade, is shown in Figure 4. 





Figure 4 


Here, the applicable equations are, 


Dor ” “29n ~ “19209 ~ “Gor Sa 
ee eo) ae toe) een COI) 
ORE blot To SiGe aR OS) 
Bey  Tl2oR "Toon: “Son eee ee) 
where Ki2 and R34 are as originally stated, but, 
Sok) TUMOUR Cae) 
Bm = Vooy (1026) 


as given by the boundary conditions. 


Thus, the finite differ- 


ence equations for the flow velocities at the bottom of a 


Solid boundary or cascade blade are, 


_ Seas 
22R 


‘ SGT OSAVAIR 


ak J+2B,[K 


L121 pare OSV 


2 2 2 
C1 A,8,) a (2B) 
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1tBy* 


(103a) 





CL-AS BL) IK, 6p tk 3,5 tBy Ky op)-2By [Ks 6p tka, p 7B Ky 07) 


U = 
Z21 2 2 (103b) 


yee 12R ~ "22r tT 87" 20% (ieee) 
Mee 121 ~ “o0r ~ “2"20R (1034) 
Voor and Voor are given from the boundary condition. 


To obtain the conditions on the initial left-running Mach 


line, assume that Poi and Poi as shown in Figure 2 are just 


in the freestream and then let Ax shrink to zero. Since 


Ki = Ke = 0 and A, and By go to zero as Ax goes to zero, 
Mies initial finite difference equations for — = constant are, 
U,> as Co5 = 0 (104a) 
U,5 - V55 A Vio = 0 (104b) 
WV = 
U,5 1 T99 0 (104c) 
Thus, 
ape Yop = & 99 Sa 


Eq. (90a) can then be written, 














(S%) = (106) 
x 
E 
Integrating gives, 
M2 
U = Ujexp[-ik x ] (107) 
- M?-1 i 
x=0 
Thus, 
uU (0) eae) (0) si He x) (1084) 
DOR 7 cos ey O)sin (k x) (108a 
R Y2oR M2- 1 IT ae 
U mu? . M? 
2 a “Vo97 60) cos (k . ; x Voop 60) sin (k ee f x) (108b) 
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Voop = ~- U 


22R 22R 
Yoo, ~ “bea 
Coon = ~ Voor 
aoo) ee Say 


22 Jaga M 


For the initial right-running Mach line, Pil 


Cireigic ) 
(108d) 
(108e) 
(108£) 


12 of 


Figure 2 are just in the freestream and Ax is brought to zero. 


Since Rio = Key = 0, and A, = BL = Q, 

UD. pee © 

sg mimelbis) <2 2 

Ujpaine yetpeae 
Thus, 

ee aa 9) uaa) 


Eq. (90b) can then be written, 

















rea = - ikU 
Integrating gives, 
M2 
ik 
u = vo) e!~4 ae a 
Thus, 
(0) cos( see (0) uM" 
U = V cos(k x)+V 0)sin(k x) 
Wi, 
2R Om oe 22M on 
M? M? 
U5 517 Vo97 (9) cos (k ; x) -Vo op (0) sin(k ; x) 
M*- M“- 1 
Voor ~ Voor 
Coon Moor 


36 


(109a) 


(109b) 


(109c) 


(110) 


Clit) 


(112) 


(113a) 


(113b) 


(113c) 


(113d) 





Voor ~ Yaar Ciise) 


Coot = “Voor CLis£) 


The pressure distribution along the surfaces of the blades 
can be determined from Eq. (42). Dividing both sides by p, 
and using Teipel's amplitude function for C, the equation 


can be written, 


(Ss ECs 59) Re at (114) 
x=0- 
where 
P wee. 
P(x,y)e = 5 


Thus the non-dimensional pressure distribution along the 
top surface and the bottom surface of the cascade blade can 
be computed in a point by point calculation of U, V, and C 


mrome lines of — = constant and n = constant, 
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IV. ELEMENTARY THEORY FOR SLOWLY OSCILLATING CASCADES 


Sauer'’s (1950) solution for an airfoil oscillating at 
low frequencies in an unbounded supersonic flow can be applied 
to the flat plate cascade to form a theory suitable for con- 
parison with the method of characteristics. As shown in Gar- 
rick (1957), the equations of motion of supersonic flow over 
a flat plate oscillating at small amplitude can be written, 
assuming irrotational flow, as 
a ae =o, 20 (115) 


Cc 
co 


cot7ad =® 
xX 


With the assumption of simple harmonic motion, 


O(x,y,t) = o(x,y,kjer** | (116) 
and the potential equation becomes, 

cot*ad = ot 2ikM’$ - k?*M*o = 0 (117) 
where 9 is complex. 

If the frequency of oscillation is Seger iene low, the 
perturbation potential amplitude may be expanded in a Taylor 
series, such that 

o(x,y,k) = x(x,y) +k p (x,y) (118) 
neglecting higher order terms. Eq. (117) then splits into 
two simultaneous equations since k cannot appear in any rela- 
tion, it being independent of k. These equations are: 


cot*a a: See (119a) 


2 _ _ a ay2 
cota aa Voy 2iM Kes (119b) 


38 





an Se 





Sauer (1950) showed that a general solution to these two 


equations is, 





y Gai CO), (120) 
WS eeZ eet a yg(Z) G21.) 
Z=x- y cot a (122) 
UD) ie reek) (123) 
Zo - So + fy Coe (124) 


where g and h are arbitrary functions of position and are 
equal to zero for z~S 0 (x,y, z and d are non-dimensional). 

To apply this general solution to the cascade with super- 
sonic leading-edge locus, the flow field is divided into 
Separate zones as shown in Figure 5. The number of these 
mees 1S dependent on the size of A. For A > 1, there ts 
one zone along each blade; for 0.5 < AS1, there are two; 
Mor 0.33 < AS 0.5, there are three; for 0.25 < AS 0.33, 
there are four; etc. 

If oscillation in pitch only is analysed, the boundary 
conditions on the blade surfaces can be written, in non- 
dimensional form, from Eq. (33) and Eq. (35), 


$, = -1 - ik(x - x) at y = 0 (125) 


o- 
i 


So Roe kee x0) ac ye 6 (126) 


where oes is one. 


39 








Figure 5 


In Zone I there is no interference from the upper blade, 


hence, Sauer's solution for the single airfoil applies. The 


boundary condition, however, must be satisfied at y = 0, 
hence, 
= = ' — 
Is cot a g (Z) 1 (127a) 
M 
= = ' SS =— fj = 
Sey cot a h (Z)-i se g (4) Cx x) (127b) 
This gives, 
g, (4) = Z tan a (128a) 


ho (2) =-i tan aZ (x + = tan*a) (128b) 


where, 
6 = g)(Z)+ kfh,(Z)-i —o— y 2) (2)} (129) 


In terms of Teipel's amplitude functions the flow properties 
can be related to those obtained in Method of Characteristics 


BY , 
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U = 9 (1304) 





x 
Y= tan @ oe (130b) 
Cee |e + ikd] (130c) 

In Zone I, 
Ul= tan a {1 - ik[x + x tan*a + merce tec | } (131a) 
y;= tan a {-1+ ik[x - x + y(cot a + — ~)}} Gem 
Cy = tan a {-1+ ik[x + x(tan*a-1)+ 2y Teor 1a) t (131c) 


In Zone I' the perturbation potential can also be written 


directly from Sauer's single airfoil solution, 


M 
cos a 





@ = ¢,(Z) + k {h,(Z) + i y g,(Z)} (132) 


where, 


Ni 


=x -A- B+ y cota (133) 
This ensures that Z < 0 upstream of the leading-edge of the 


upper blade. 











The boundary condition at y = d gives, 
= Tso es ié 
Xy cot a g, (4) e (13462)) 
= 1(7 ‘ 7 ° 1(7 
eS cot a h (2) + eras g(4)+ Lae scot agi (2) (134b) 
hence, = - i(x - Xo 7 B) eis 
g, (2) = - Z tana ane Ci 352) 
eee - id Z 2 Md 
h (4) = itan a Ze Wes) 5 tan“at Tears (135b) 
Then, 
Uy -tan a e*°{1-akIx +(x-B) tana +A-y cot aj} (136a) 
Ss ore ah! _ oM 
Vie Ban a e  { 1+ik[x, (x-B)+(d-y) Ccotat ae ada (136b) 
_ -tan Qa eo °{-1tik x +(x-B) (tan?a-1)+2A-2y ect al} (€i36c> 
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i. 0.5 = A = 1.0, Zones Il and Til" will be present in 
addition to Zones I and I’. In Zone [II the perturbation 
potential is due to the initial waves from Zones I and IL' 
mes the reflected waves at Zone ITI, hence, 

_ g (4) + g, (4) a g, (Z,)tkih, (Z) th, (2) 


VAeeeZ = gmeDiaHe, (4,0) (137) 





; M 
. h,(Z,)-2 cos a 


where Z and Z are defined in Eq. (122) and Eq. (133), and, 


Z4 =x - 2A + y cot a (138) 


The boundary condition at y = d gives, 


SS 7 id 
Ky cot a g.(2)+ cot a go (2) + cot a g, (2Z,) =e (139a) 


=—_- t t 1/7 ae 
ie cot a hi (Z)+ cot a h i (Z)+ eoe o h, (Z,) 


M 
cos a 


Md 
cos @ 





- i le (Qe Bo Caos 3a He) 





Co tne [gi (Z)+ BC) se g,(Z,)] 


id 
a (x -(xXe -Eee (139b) 


hence, after some algebraic manipulation 





g, (Z,) : Zz, tan a (140a) 
a ; 1 Md 
h,(Z,) = -i tan a Z,ix + 7 tan*a + 2 ae au (140b) 
Then, 
See— tan ai2 - eid -~ik[ (2 i6 Z 
II i -e ) (x + x tan“a + A) 
id 2 
jee (yecot Oo +7 Betanco me (141a) 
V..= tan a mero, Lk i6 - - a 
Il { iteie (x, x +(d-y) (cot at ae OD 
-2(d-y) (cot a + ——) ]} (141b) 
Cio -tan Mims Ok ikfe*® (x +(x-B) (tan?a-1) | 
+ 2A - 2y cot a)- 2(x + x(tan*a-1) + 2A)]} (141c) 
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In Zone II' the perturbation potential is due to the ini- 
tial waves from Zones I and I' plus the reflected waves at 
Zones II and II', hence, 
pes, (2) + 2, (Z)+ 8, (Z,)4 kth (Z)+ h (4) 


pez) - 4 yle((2)= 2,@)+ 2,@,)13 (142) 


ee 
cos 


where Z and Z are defined in Eq. (122) and Eq. (133), and, 














Zs = (x - A - B)- y cot a (143) 
The boundary condition at y = 0 gives, 
Bay ’ ‘re 
Xy = cot a g (2) + cot a g, (4) 
== 1(7 = = 
cot a g; (2,) 1 (144a) 
a ! 1(7\— 1/7 
Vy cot a hi) (2) + cot a hi (Z) cot a h, (2,) 
-i ——"— [g (Z)- 2, (Z)+ 2,(2Z,)) = i(x ~ x,) (144d) 
hence, 
sia, = Shzelwwaneme° Gace 
fe Bg! 
Z, 
ae =| id ee 2 Md 
h, (2,)= itan a Z e {x0 + 5— tan“o ae ay (145b) 
Then, 
7 id : i6 2 
Uyp.= tan a{1-2e + ik[2e (x +(x-B) tan OL 
+ A) - (x + x tan*a + y cot a)]j} (146a) 
e sia id M 
tee tan af-e “+ ik[2e”(y(cot a + nae =e, 
M 
- (x o> re he S@evaie Ch 47 re ee (146b) 
Cry -tan a{ Lome ik[2e"° (x +(x-B) (tan?a-1) 
eee Cat axcane ce) ey cor (aya (146c) 
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Because of the algebraic complexity of further analysis, 
additional zone calculations will be made with the assumption 
that the blades are oscillating in phase, i.e., 6 = QO. 

Hote. ss = A = 0.5 in addition to Zones I, I', II, and 
II’, Zones III and III' will be present in the flow field. 

In Zone III the perturbation potential is due to the initial 
waves from Zones I and I' plus the reflected waves at Zones 
mee Ii’, and III, 

@ = 8,(Z)+ 8,(Z)+ 8,(2,)+ 8, (2,)+ 8,(2,)+ 


k{h(Z)+ ho (2)+ me AS cr h,(2,)+ h,(Z,) 


meaty _ yee : 
1 tos a8 (4)- 8) 0Z)- 8, (2,)+ 8, (2,)+ 8,(24)]3 (147) 

where Z, Ze Zi» Z4 are defined as before, and, 
Z. =x - 2A - y cota (148) 


The boundary condition at y = 0 gives, after considerable 


algebraic manipulation, 





85 = Zo tan a (149a) 
> ‘ Md | 
he = - tana Zo {x + 7 tan°a + 2 aos a! (149b) 
Then, 
Vir = tan a{l- ik[x +(x-2B)tan*a + 2A+ y cot aj} C1503) 
a8 dae - M 
Vig =-tan a{l ik[x, x + y(cot a + aoe ie (150b) 
Cror =-tan a{l- ik[x +(x-2B) (tan*a-1)+4A t2y coma), CuS0c) 


In Zone III' the perturbation potential is due to initial 
Waves from Zones I and I' plus the reflected waves at Zones 


Mie ll’, and IPL’, 
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= 8 C4)+ 8, (Z)+ 8, (Z,)+ 8, @,)t 8,4) 
+ kth (Z)+ ho(Z)+ hy (Z,)+ by + bL@,) 
my oy e 7) - 7) 7 
Le ol 8o 62) - 8, C4)- 8, C2404 8,04) g,(Z,)]} (151) 
Piere Z, Z, Za» and Z4 are defined as before, and 
Zo =x - B- 3A + y cot a Gs) 
The boundary condition at y = d gives, after considerable 


algebraic manipulation, 








SpeS ~ Ap bemac (153a) 
q rs my) ; Md 
h, (Z,) = itana Zo {x + 7 tan a + 3 ae a C1535) 
Then, 
Urq - tan afl - ik [x +(x-3B) tan* a+ 3A - y cot aj} (154a) 
; M 
Vari - tan a{l-ik[x -(x-B)+(d-y) (cot Qa + Fea me (154b) 
Cyyr'" tan a{1-ik[x +(x-3B) (tan*a-1)+6A - 2y cot al} (154c) 


more 0.25 => A S OF33emin addition to Zones 1, 1I', Tiare = 
IIL, and III", Zones IV and IV' will be present in the flow 
field. In Zone IV the perturbation potential is due to the 
initial waves from Zones I and I' plus the reflected waves 
Memzenes II, IL', III, II1', 
$= gi(Z)t+ B(Z)+ 8, (Z,)t+ 8, (Z,)4 8, (Z,)4 8,(Z,)+ 8, (25) 
+5 k{h (2) + h (2) + hy (20 F hy @,)+ hy (Z5)+ h, (Z,)+ h,(Z,) 
a a oP) 7 m7) 2 
45 Ig (2)-g, (2) -g, (2) +84 (Z, +8, (24) By (Z,)-8 4 (25) )} 
(1855 ) 


where Z through Zo are defined as before, and 


Z3 =x - 4A + y cot a (156) 
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The boundary condition at y = d gives after considerable 


algebraic manipulation, 


_ 8, (23) = Z., tan a 





a 2 Md 
h,(Z,) = -i tan a Z,{x + 5 tana + 4 re a 


Then, 


Uy tan af1-ik[x +(x+3B)tan*a + 3A + y cot al} 


M 
Cosmd 


)}} 





Vive -tan a{l-ik[x -(x-B)+(y-d) (cot A+ 


Crys -tan afl-ik[x +(x+3B) (tan*a-1)+6At2y cot al} 


In Zone IV' the perturbation potential is due 
tial waves from Zones I and I' plus the reflected 
momes Ll, ILl', LIL, III‘, and IV': 
¢ = gi (Z)+ Bo (Z)t+ 8, (Z,)t 8, (4,)+ 85 6Z5)+ 8, (25)4 

tkt{h (2Z)+ h (Z)+ hy (2,)+ h, (2,)+ h,(Z,)+ h,(Z,)+ 


My 
cos a 





[g,(2)-8,(2)-8, (2, ) +8, (2, ) +8, (Zy) 
-g,(Z,)+ @,(Z,5)1} 


Where Z through Z, are defined as before, and, 


2 
Z. =x - B - 3A - y cot a 


The boundary condition at y = 0 gives, 


g,(Z,) = Z., tan a 


Ni 


Md } 
cos a 





= x: = me) 2 
h,(Z.) aeomeanl G2, {x + tan‘a + 3 


KO 


Then, after considerable algebraic manipulation, 


Uryr= - tan o{l-ik[x +(x-4B) tan*o + 4A - y cot al} 


—M 
cos Qa 





) 3} 


Viyie - tan at1-ik[x)- x - y(cot a t+ 
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C57 a) 


(157b) 


(158a) 


(15 8b) 


(158c) 


to the ini- 


waves at 


g,(Z,) 


Ma (2a) 


(159) 


(160) 


(161a) 


(161b) 


(162a) 


(162b) 





C tan a{l-ik[x -(x-4B) (tan*a-1)+ SA = -Zy cot ol) Cle Ze) 


a 

Further zones may be calculated for values of A less than 
0.25 in the same manner. 

In the Method of Characteristics it was shown that the 

non-dimensional pressure, P, is a direct function of C. 
By integrating this pressure over the blade surface, the lift 
and pitching moment on the blade may be determined. The deri- 
vation of the lift and moment equation is given in the Flutter 
Analysis. In terms of Garrick and Rubinow's (1946) lift and 
moment coefficients, for the case of in-phase blade oscilla- 


tion (6 = 0), 


9 1 a: i+B a 
L=L, +itL, = ~- —{f C(x,0 )dx - fs C(x,d )dx} (163) 
3 4 2 
k* oO B 
: 4 : + 
M=M, + iM, = - —{f (x-x ) C(x,0 )dx 
3 4 y. O 
k? © 
1+B - 
-~ f (x - B- x0) C(x,d \dx} (164) 
a 
For 0.25S AS 0.33, 
i 70 
L, = ge ano td2 ereoke 4 Bo (tan coon , Ces 
M, = -16 2 gap el AON lag (166) 
k 2 
ay 16 - tan a{x (d?+2a7+B*tan*a)-8A°~4d7A+4AB?} (167) 


If the cascade blades are allowed to oscillate at an ar- 
bitrary interblade phase angle, the pitching moment determi- 
nation can be somewhat simplified by integrating in the 


Following manner: For 0.55 As 1.0, 
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“a A-B = 
{ i Cx-x J 10, G,0 ) 





16 


- 3S A 
=e C,1(xtB,d YJdx + f 


+B ne 
A-B 


7 aid C__(x+B,d )Jdx + , es ue (x oF) 
Il ; A+B ns 


Ba? C1, (x+B,d )]dx} Clore) 


This integration gives 


M, = 8 tan afcos 6 [x (2-2A)+(A-B) *-1-2AB] -k B sin 6 
k 
ex: a eee cade) WA — (range ) 
(A74+ 1 - J R®y } - x + =a (169a) 
3 re) 2 
p, = 8 = tan a{cos S [x? (2A=2)+ CG — 4A + 2 


Me 


fc’ = B°- L)tanea)- =~ A*+ 2¢1-B°)A = 5 


B sin 6 [x_ - A] 


2 (1-43 Qe) vs 
+ 3 (1-A°)tan‘“a] 2 m 


a 
Ms 
N 
[ 
Ms 
+ 
| 
a 
oo 
| 
s 
J 


=) tan*a} | (169b) 


In most airfoil theory, the non-dimensional pitching 
moment is written as, 
C =68[c + ikc .Je*t (170) 
m o- MM, m6 


Relating this to Garrick and Rubinow's (1946) moment coeffi- 


cients gives (where 0 is one), 


Cc = - Sk? M (171la) 


(171b) 
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VY. FOURPERP er yrels 


Lane (1956) showed that the flutter analysis of a cascade 
can be made by considering only a single blade of the cascade 
oscillating in the flow field. In order to compute the non- 
dimensional lift and moment acting on that blade, the pressure 
difference over the top and bottom surfaces must be integrated 


over the blade. Thus, 


L = fe (ees 0u) ~ P(x,0 )}dx (172) 
O 
1 
M= ff (x- x) {P(x,0")- P(x,0 )}dx Clg ) 
O 


where P is the mOnmdeine née oma pressure as defined by Eq. 
rid). 

Lift here is positive downward as in Garrick and Rubinow 
(1946). The non-dimensional pressure on the lower surface of 
this equivalent blade is that of the lower surface of the 
blade above adjusted for the interblade phase angle and 
stagger such that, 

De Oe ES ae P(x+B,d) enZ> 


The lift and moment on the equivalent blade can also be 


written 
: + at ee 
ieee r(x, 0 )dxea a) e P(x,d)dx (175) 
. B 
; + Zine -i6 
Me= f (x-x 9) P Cx, 0 ) dx a (x-B-x ))e PCa) dx 


i (176) 


separating these into their real and imaginary parts and writ- 


ing the equations in terms of the local sonic velocities gives, 


49 





1+B 
{[cos Q°C, (x, d) + sind*C,(x,d) dx} 
(177a) 


1 
a 2 + _ 
Pe ye if C,(x,0°) dx ‘ 


Bee Rn 1+8 
Ly= YM {f C(x, 0 )dx - f[cos §°C,(*,d)-sin 6 °C, (x, d) dx} 
(0) 


: (177b) 
. 1 ae 1+B5 
ee YM if (x-x) C, (x, 0 ) dx a (x-B-x ) 
[cos 6 C,(x,d)+ sin 56°C, (x,d) Jdx} (178a) 
A 1 a. 1+B 
= YM ud (x-x)C,(x,0 )dx a (x-B-x ) 
[cos 6°C,(x,d)- sin 56°C, (x,d) ]dx} (178b) 


where Ce and Cy are the real and imaginary components, res- 
pectively. 
To evaluate these integrals numerically, the trapezoidal 


rule is used: 


a £f(x)dx = Sx LE (x) + 2£(x), + 2£ (x) 5+ 2£ (x), ae 


re ae = 2£(x)_ot 2£(x)_4+ f(x) ] (179) 


If oscillation only in bending and torsion (pitching and 
plunging) is considered and Garrick and Rubinow's (1946) 


method of expressing lift and pitching moment is used, 


2h 
me. Ll A 90, 6 O ‘ ° 
16 7 Po © Uy k“e ie eis ig aL L,)t+ 8 (Lat il L,)] (180) 
Ms - i eT pPeWt Oy Gade Ma) ORR ciate) Cla) 
4 Peo co a es eZ Neg itary 
Where k = unlike Garrick-Rubinow's, ts mothe yertical 
Lee) 


equilibrium position of the elastic axis, and oe is the 
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equilibrium angle of attack of the blade. In dimensional 


form the lift and pitching moment equations are, 


L=@p,e  {[L,+ £ by] +IL,+ £ Ly] } (182) 
Plunge Pitch 
mode mode 
& 
fe= e?p ec * {{M,+ i M,] Altera Mend } (183) 
Plunge PAE Cit 
mode mode 
Thus, 
Cc 
L, t iL, = - —— [L, +i L.] (184) 
. : M*k7y h : Plunge 


Separating the real fr 


Nn 


Cc 











mode 


om the imaginary parts of the equation, 


a im a 21.2 Le Plunge (185) 
M“k“y h 
oO mode 
OY, 
L,= ~ {fi -Cx(,0 yax — fl eos 6*C7G.d)tsin 6 
2 R R 
ke he B 
'e) 
Cy (x,d) Idxt a vice mode (186a) 
Likewise, 
a 1 fe 1+B 
L,= - {i SCRE idx — Jl COSHOURIMICEGera) -— sinec 
eo ere : B ‘ 
’e) 
"CR (x,d) Idx} 51 un ge mode see 2D) 
1 1 2 Raw 
L,= -2 io EGO) dso eeros.6 -° mG (x,d)+ sin 6 
3 2 R R 
k* 6 Oo B 
’e) 
Cy Ce, a) dxd os ay mode ES) 


Sal 




















1 J + spt 
a =? {f/f C1 (x, 0 \dx - { [cos §+C) (x, 4) 
[eo a, 9 -B 
- sin 6*C, (x, d) J dx} (187b) 
Pitch mode 
a 1 a: 1+B 
b> 2 {f (x-x )C, (x, 0 ydx - f (x-B-x_) 
ko hf oO x B : 
Oo 
*[cos 6 °C, (x,d) + sin 6*C,(x,d) Jdx} 
Plunge mode (188a) 
a i+B 
> -—2 tr (x-x Pie, (x, oF ydx - ae (x-B-x -, 
k* ho 
O 
[cos 6° “Cy Cor sano eG pbXs d) Jdx} 
Plunge mode (188b) 
1 1+B 
M,= —4 tf (x-x 5 SR (x, oF ydx - f (x-B-x >? 
k? @ Oo. B 
Oo | 
[cos 6°C p bX: d)+ sin 6° Ci (x, ance: 
Pitch mode (189a) 
| 1+B 
M)= -4 ir (x-x ir. (x, ot ydx - f (x-B- = ) 
ke cick ° B 


*[cos 6 °C, (x,d)- sin 6°¢C Pes d)}jdx} 
Pitch mode (189b) 


where the integrals are computed using the trapezoidal rule 
and the two different oscillatory modes are determined by 
the boundary conditions at the blade surfaces. 

Bisplinghoff, Ashley, and Halfman (1955) define flutter 
as "the dynamic instability of an elastic body in an air- 
Stream." In Chapter 9 they describe the physical nature of 
the flutter phenomenon. Basically, flutter occurs when the 
piaeveristic determinant of the equations of motion of the 


blade vanishes. 
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If blade oscillation in two degrees of freedom (pitch and 
piunge) is considered, the equations of motion of the cascade 


blade can be written as, 


‘Mh + SV t Ch = L (190) 

S gh a 1,8 + CA = My (191) 
where the total inertia force acting on the blade is, 

Fy = - (Mh + S99) (192) 
the elastic restoring force is, 

Pe = = Cyn (193) 


the sum of the moments of the inertia forces about the 
elastic axis is, 


My = - (S gh =) I,8) (194) 


and the elastic restoring moment is, 


M. = - C,6 : (195) 


L and My are the aerodynamic lift and pitching moment about 
the elastic axis. With the assumption of simple harmonic 


oscillation, the blade motion can be written, 


Lot 
e 


aS. (196a) 
h = - w? we (196b) 
8 = 6. ae (197a) 
8 = - w? 6) ais (197b) 
The equations of motion thus become, 
= LWt a 2 _ 2 
iE = e [- w°Mh, - w S49 + cho d (198) 
M.= 2 w*S.h - wT Ce tr eGcael (199) 
Q Go 8° o 6-0 
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Equating these to Garrick-Rubinow's standard lift and moment 
expressions and following Garrick-Rubinow's (1946) method, 


the equations can be written as, 





2h 
X=— es = = 
(2) ue L,t L L,) 5 + ( UX gt Lt 3 ee 0 (200a) 
2h. ; 
(-UX, + M + i a + (2X = UX 5 + M, + i M,)80= 0 (200b) 
The two simultaneous flutter equations form an eigenvalue 
2h 
problem in and OS? thus, the determinant of coefficients 
Cc 


(the characteristic determinant) must be zero: 


2X - pth + 2G — UX, FL i LG 


al Ze 3 4 


= 0 (201) 
tM, + i M, 2 4X- ure oe Jal ae hy 

This is the flutter determinant with unknown X. If the 
determinant is expanded, and the real and imaginary parts 
equated to zero the result is two equations that given a 
value of k and M will yield values of X. For a specific 
Mach number, the flutter point occurs at the value of k that 
produces identical values of X in the two equations. From 
this value of k and X the flutter speed and frequency can be 
calculated. 

The two flutter equations are: 
(Real part) 
MM pXPF[My(Ly-w) + A Chz-urg)]x + Cy = 0 (202) 


(Imaginary part) 


(QqL. + QM, )X + Cy = 0 (203) 


where, 
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S meee 2 oe 2 
= = = 2 


and where, 


D a => Lele - Lam 


Be ys seed ee 


uM, (205a) 


D 


I LM, = LM, e LoM, = LM, 


(205b) 
Solving for X, 
(Real Part) 


-[2, (L,-H)+2, (M,-yr 9) 14/72, (Ly -n) +2, (nur, 1 °-4C, 0, 2, 


- : R“h''6 
ne (206) 
(Imaginary Part) 
7 GY 
eC) (207) 
24 Lo + 2) M, 


The solution to the flutter problem is then to obtain 
values of x as a function of k until the real and imaginary 
solutions converge. At this "flutter point’ the non-dimen- 


sional flutter frequency is, 


ff 


yee 


= (208) 


and the non-dimensional freestream flutter speed is 


ea 


Q 


(209) 


Q> 
: 
=I 


BREvshculd be noted» Ehat@inestnte blade spuase angle a son 
open parameter (as well as k) in calculating lift and moment 
on the Bae cale blade. Lane (1956) showed that the minimum 
flutter speed of a cascade occurs at a non-zero interblade 


phase angle and that this is the critical phase angle at. 
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which the cascade will oscillate. For an oeaeecved cascade 
this. critical phase angle was shown to be 180°, but for a 
staggered cascade it can only be found by repeated flutter 
calculations at various phase angles. The problem is solved 
when the minimum flutter speed is determined. 

If only a single degree of freedom of motion (pitch os- 
cillation) exists, the flutter problem is greatly simplified. 
The equation of motion of the blade can then be written, 


I,6 + (1+ig)C,0 = M, (210) 


Here the sum of the moments of the inertia forces about the 
elastic axis is, 


M, = = 1,8 (231) 


the elastic restoring moment is, 


M, = - (1t+ig)C,8 (212) 


and the aerodynamic pitching moment about the elastic axis 


iS, 


a: i 
Hy SR he 


A2 525.2 Wt 
soe is Oe (M+ at M,) (293) 


With the assumption of simple harmonic oscillation as before, 


the blade motion can be written as, 
@= 6 et (214) 
— é ene (215) 
The equation of motion then becomes, 
eNO [-Tgu? Ot (itis) C0, )= - Fhe u2k7O 2 ME Cti My) (216) 


Separating the equation into its real and imaginary parts 


gives, 
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i 


Sapa: rec, + GPaoo uM, = 0 (217a) 
g Cy + Gat w*M, = Q (217b) 
or as in Garrick and Rubinow (1946) 
aX = Ur, + M, = 0 (218a) 
M, + g MX = 0 (2180) 


Thus, the solution to the single-degree-of-freedom prob- 
lem is to obtain values of M, as a function of k until the 
Second equation is satisfied. The value of M. for that k is 


then used in the first equation to determine if X is reason- 


able: 


=. 16) (219) 


os (220) 


ele 
<i 
rs 


es ee (221) 
Gd) 


Again it must be noted that in the case of staggered 
cascades the flutter computations must be repeated with vary- 
ing interblade phase angle until the minimum flutter speed 
is found. This will give the critical flutter condition. 

As in the analysis of the single airfoil, in the absence 
of structural damping, aerodynamic instability (flutter) 


occurs when = 0. When M) is written as the pitch damping 


a7 





derivative, Eq. (171b), it should be noted that aerodynamic 
undamping occurs when the pitch damping derivative is posi- 


tive. 
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~VI. COMPUTATIONAL PROCEDURE 


The computer programs used to calculate the cascade flow 
field are based on one used by Platzer and Pierce (1970) to 
calculate pressure distributions on an airfoil oscillating 
with wind tunnel interference. The programs here, however, 
in addition to calculating the pressure distribution on the 
blade of a cascade integrate the pressure distribution over 
the blade surface and calculate the flutter frequency and 
freestream flutter velocity. One program (Program A) is for 
two-degree-of-freedom (pitch and plunge) flutter calculations 
while the second (Program B), merely a simplification of the 
first, is for single-degree-of-freedom (pitch) flutter cal- 
culations. 

Both programs are valid only for cascades with supersonic 
_sadsinioedleg (Oca Vou reine Goel smell), amaieeedle egal les 
tions. The programs calculate the flow field between two 
adjacent blades using the method of characteristics finite 
difference equations developed previously and integrate the 
pressure distributions on the blade surfaces using the trape- 
zoidal rule. 

The computational molecules shown in Figures 2, 3, and 4 
are used to calculate the flow field quantities at Poo akg 
the general flow field, on the lower blade, and on the upper 
blade, respectively. The distance Ax shown in the figures is 


determined by, 
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Ax = (222) 


2A 
as 
Where v is the grid fineness ratio. This parameter is an 
arbitrary input variable into the programs and is equal to 

one less than the total number of grid points on any Mach 

line that runs from one blade to the other. 

Input parameters are entered into the programs by means 
of the NAMELIST option of FORTRAN. In using this option, in- 
put parameters are combined in a titled list and then referred 
to in the programs by this title. In both Program A and B 
the NAMELIST title is the same (NAM1); however, the input 
variables in each are slightly different. The real advantage 
of using NAMELIST is that each data card need only contain 
the value of one input variable while the format of the data 
card is such that the variable's name is punched on the card 
as well as its value. In this way, when a different value for 
an input variable is desired, rather than changing the entire 
data deck, only one easily identifiable data card need be 
changed. 

The format for a NAMELIST data deck is as follows: The 
first column in each card of the deck is left blank. The 
first card starts in the second column with the symbol & 
(ampersand) followed immediately by the title of the NAMELIST 
(here, NAM1). The next cards list the variables and their 
input values, one per card, in any order. The format here, 
Starting in the second column and followed by a comma, is: 


Variable Name = Value. The last card in the data deck and 
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the signal to the computer that the NAMELIST has ended again 
starts in the second column with the symbol & (ampersand) 
followed immediately by END. 

The value of an input variable may be written in any for- 
mat as long as integer variables are written without decimal 
points and real variables are written. with a decimal point. 
Blanks are taken as zeroes, but a comma must appear somewhere 
between a desired input value and the next data card. 

In both programs, prior to the NAMELIST dataset, the date 
is entered on a data card in the first twelve spaces. 

To illustrate the input procedures, the following two 


examples are proposed: 


awo Degree of Freedom Flutter One Degree of Freedom Fiutter 


v = 100 v = 100 
ie yo 
y=1.4 y=1.4 
k = 2 k = 2 

Xo = 0.5 Xo = 0.5 


QA. 
Il 
© 
oO 
Qu. 
tt 
© 
ON 


B = 26.565° B = 26.565" 
6 = 180° § = 180° 
oS FOC yu = 500 

a = 0.5 Xo = 0.5 

Xe = O.L g = 22 

oh 

— = 1.5 | ee 

0 

A2/k = 1 Coen 

€ = 1.2 
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The data cards needed are as follows: 


Card — Program 


ul 


2 


10 
dak 
2 


13 


13 


14 
15 


16 


16 


17 


A 


A 


> PF PY 


> 


& 


& 


B 


B 


Format 


6 WUE 197 2 


&NAM1 

NGRDFN = 100, 
FSTRMN = 1.414214, 
RTOSPH = 1.4, 
REDFRQ = 2.0, 


XSUBO = 0.5, 


TNWDST 


tt 
© 
ON 
we 


STGANG 


tt 
No 
ON 
Ww 
oN 
Wi 
ee 


FAZE = 180.0, 
MUU = 500.0, 


RSUBA = 0.5, 


XSUBA = 0.1, 

GSUBA = 0.02, 
INCRE = 1.0, 

FIN = 1.2, 


HAFREQ = 1.5, 


& END 


& END 
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“VesiecEiptlon 
Date 
NAMELIST title 
Grid fineness ratio, must 
be an even integer, less 
than 400. 
Freestream Mach number 
Ratio of specific heats 
Reduced frequency 
Elastic axis position 


Blade distance 


Stagger angle, tan BS 
COE: 


Interblade phase angle 
Blade density parameter 
Radius of are 
Distance from the blade 
elastic axis to its center 


Of Je ray ity. 


Structural damping coetiae 
cient 


2/k increment 
Chord lenerh. c= 2A, 


Ratio of natural frequencies, 
bending to torsion. 


End of NAMELIST 


End of NAMELIST 


a 





The definitions of these parameters are giyen again in the 
alphabetical listing of all program yariables in Appendix A. 
Program A has eight subroutines. They are as follows: 


INPUT Reads in all input data and sets up the 
finite difference grid. 


INTIAL Initializes all the flow quantities at the 
leading-edge, starts the integration, and 
initializes most of the logic variables. 

MACHLN Computes the flow quantities along the 
initial left and right-running Mach lines 
at the given grid point. 

HIFOIL Computes the flow quantities on the lower 
surface of the upper airfoil at the given 
grid point, and continues the integration 


of the pressure there. 


GENFPT Computes the flow quantities at a general 
flow Field erid point. 


LOFOIL Computes the flow quantities on the upper 
surface of the lower airfoil at the given 
grid point, and continues the integration 


of the pressure there. 


COMPXY Computes the x and y position of the given 
Sirid peinme. 


FLUTER Computes the oscillation frequency ratio, 
X, for both real and imaginary flutter 
equations. 
Program B contains all of these subroutines except FLUTER. 
Since the flutter problem is greatly simplified in the single 
degree of freedom case, the flutter calculations are made in 
the main program. 
Program A contains at least two "DO" loops. The inner 
pene is for calculating the flow field watiecarst pitch coca 


lation boundary conditions and then plunge oscillation 


boundary conditions. This loop is not present in Program B 
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Eemce oOmly the flow field due to pitch oscillation is desired. 
The next outer loop is used to increment 2/k (l/k in Garrick- 
Rubinow's (1946) notation) starting at that reduced frequency 
specified in INPUT. This loop is continued until hopefully 
the flutter point is found. Any loops outside of these are 
used to make multiple flutter calculations in one run with 
Varying input parameters. 

In INPUT the grid is determined using the grid fineness 
ratio, blade spacing, and Mach number. The value of Ax as 
Shown in Eqs. (96) is obtained. The stagger angle is then 
adjusted such that the distance the upper blade is staggered 
back is the nearest multiple number of Ax increments. The 
"Compatible Stagger Angle" is then printed out. The upper 
blade may be staggered back such that the distance between 
its leading-edge and its intersection with the initial left- 
running Mach line from the lower blade is zero, but no further. 
With less stagger this distance must be greater than or equal 
pe 2\x, or, 

A Bee | (223) 
This allows a minimum of three grid points in the primed flow 
field zones as shown in Figure 5. This is the minimum number 
of points required when the value of the last point must be 
linearly extrapolated from the other two, as is the procedure 
pm the programs. 

Flow quantities along the initial left and right-running 
Mach lines are computed using the equations derived in the 


Method of Characteristics for initial conditions. Those 


64 





quantities along the blade surfaces are computed using the 
equations derived for the appropriate boundary condition. 
The values of the normal yelocities are derived as a result 
of the analysis of the flow tangency condition and have been 
given previously. 

Because the programs compute the flow field in rows of 
right-running Mach lines starting with points on the initial 
left-running Mach line, provision must be made for the com- 
putation of points in the area referred to as Zone I’. The 
reason is that once the computation reaches the intersection 
of the initial left-running Mach line from the lower blade 
and the initial right-running Mach line from the upper blade, 
the points on the initial left-running Mach line become 
general flow field points requiring those equations for the 
flow quantity computation. To use those equations the flow 
quantities along the last left-running Mach line in Zone I’ 
must be known. They are obtained in two ways depending on 
the value of the interblade phase angle: If the phase angle 
is non-zero, the entire flow field in Zone I' (along with the 
integration over the bottom surface of the top blade) is com- 
puted concurrently with the computation of the flow field in 
Zone I. If the phase angle is zero, only that row in Zone [I 
that is the image of the last left-running Mach line in Zone 
I' must be determined since the flow quantities of these two 
rows are the same (except for the signs of U and C). 

The pressure distributions on the surfaces of the blades 


are integrated point by point as they are being calculate 
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since the pressure coefficient at the last point in each zone 
of influence is unknown, its yalue must be linearly extrapo- 
lated from the preceding points. The pressure distribution 
of the upper blade is adjusted so as to represent the pressure 
distribution on the bottom surface of the lower blade; 

C(x,0 ) = eid C(x+B,d) (224) 
In this way the result of the integration is the lift and 
pitching moment on a representative blade of the cascade. 

Termination of the flow field calculation starts when the 
distance from the leading-edge of the lower blade to the 
point of calculation on the lower blade reaches the prescribed 
chord length. Because HIFOIL must be called at least once for 
proper program termination, this chord length has a lower 
bound : 

Buses DIA | = (C225) 
Since the flow field is being calculated along rows of right- 
running Mach lines, the last zone along the upper blade is 
then computed to complete the calculation. Although the flow 
field is terminated when the calculation reaches the prescribed 
chord length, the integration of the pressure distributions 
is terminated when the chord length reaches unity. 

When the values of total lift and pitching moment haye 
been determined, Program A enters them in the flutter equa- 
tions, Eq. (206) and Eq. (207) via subroutine FLUTTER. Program 
B, computing only the pitching moment, enters these values in 
pqs. (218). Cent cence to the flutter point is then deter- 


Mined. 


66 





in Program-«A, the difference between solutions to the 
meal flutter equation and that to the imaginary flutter equa- 
eon is noted after each computation. If the absolute yalue 
of the difference is less than Tomes the flutter point is 
considered to be determined. When the difference changes sign, 
indicating convergence somewhere in the frequency interval, 
the 2/k increment is divided by five and the computation is 
Started again at a new value of 2/k determined by adding the 
new increment to the previous lower value. Each time the 
flutter point is passed the increment is made smaller. The 
flutter point is considered to be ascertained when the incre- 
ment is less than or equal to Tome. 

In Program B, the value of X is determined from Eq. (218a). 
This value is then used along with M, to determine the aioe 
of Eq. (218b). If the absolute value of Eq. (218b) is less 
then Oe: the flutter point is considered to be determined. 

If Eq. (218b) is negative, indicating convergence somewhere 
in the frequency interval, the same procedure is followed as 
in Program A. When the flutter point is determined, the 
value of X is tested to insure a realistic flutter frequency 
as Shown in Eq. (219). 

The output from Program A consists of the input parameters 
(including such items as x and the compatible stagger angle); 
the yalues for Xp and Xy as defined Jue tawec2 uo) -and Eq. (207). 
respectively, for each reduced frequency calculated; and the 


non-dimensional flutter frequency and flutter speed as defined 


misikq. (208) and Eq. (209), respectively. The yalue of 
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"1/k" given in the output corresponds to the reciprocal of 
the reduced frequency as defined by Garrick and Rubinow 
(1946). This is equivalent to our 2/k. 

The output from Program B is identical to that of Pro- 
gram A except that Program B gives the values of My» M,; A, 
and the system damping (M, ~ gt, X) as shown in Eq. (218b) 
for each reduced frequency calculated. 


A flow diagram of Program A is given in Appendix B. 
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A. COMPARISON OF ELEMENTARY THEORY WITH MILES 

Miles Mee yiecesented an analysis of supersonic flow 
past an oscillating airfoil subjected to wall interference 
in a solid-wall wind tunnel. Using Laplace transform tech- 
niques he derived the following results for the moment co- 


Meticient (his Eq. (7.18)) 


“mo= 4 tan at (2N+1) [x - SI-N(N+1)Ax + ZN (N+1) (2N+1) A?}(226a) 


E . ae 1, 
Bc 4 tan af (2N+1) [x* x + =] N(N+1)A 
2 a rer oe } 
pcr xt al + Fon (N+1)“A*} (226b) 
c .= 4 tan eon) tance at x eNO 
ma 2 ~o 3 
[x =1+ + N(N+1) (2N4+1) (tan2a +2) 
° Nos TN? (NH1) 2 (2 tanza + 3)A°*} (226c) 


where N is the largest integer smaller than at and a, @ 
subscripts represent angle of attack. Here, for a slowly 
Secillating airfoil, 


LKE 


C = a teat I + c .)je (227) 


m ma 
unlike Eq. (170). 

As stated by Miles, this case is identical to that of 
the unstaggered cascade oscillating with an interblade phase 
angle of 180°, and the yalues of pitching moment as defined 
in Eq. (170) and Eq. (171) based on our elementary theory 
for B= 0, 6 = 180, and 0.55 AS 1.0, prove to be the same 


as those shown in Eqs. (226) with N = 1. The torsional 
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stability boundaries Sie = 0 for an unstaggered cascade os- 
cillating with an interblade phase angle of 180° as computed 


from Eq. (171b) are shown in Figure 6 for yarious values of 


solidity. ; 
A= 0.7 
Mo 
AN ESOS 
= 
1) A= 0.9 
| 
Unstable (A < 1) 
| 
D 
= 1 @ 
Unstable 
1 


afl ‘2 3 wy ‘5 





Figure 6 


The large destabilizing influence of interference causing 
instability over the whole supersonic Mach number range for 
a wide range of elastic axis positions can be seen. Also 
shown is the interference - free boundary (A = 1) for the 
single airfoil in an unbounded supersonic flow (Garrick and 
Rubinow 1946) exhibiting instability as is well known only 


for low supersonic Mach numbers. As demonstrated above, 
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these stability boundaries also apply to the slowly oscililat- 


mig airfoil subjected to solid-wall wind tunnel interference. 


B. COMPARISON OF ELEMENTARY THEORY WITH LANE 

Lane (1957) extended Miles’ Laplace transform solution to 
the analysis of oscillating cascades with supersonic leading- 
edge locus. By expanding the blade pressure distributions 
(his Eq. (13)) with respect to frequency and retaining only 
terms up to the first power in frequency, the pressure jump 
across a slowly oscillating blade can be obtained. This 
pressure jump is identical to that being integrated in Eq. 


(168). 


C. COMPARISON OF ELEMENTARY THEORY WITH DRAKE 

Drake (1956) using Laplace transform techniques presented 
an analysis of supersonic flow past an oscillating airfoil in 
a wind tunnel with free jet boundaries. His pitching moment 
coefficient is defined in the English system, equivalent to 
Ba. (170), is, 


KE 


_ a oll 
Cs = am, + ik mn. e (228) 


for a slowly oscillating airfoil where a@ again refers to 


angle of attack. Here, 


my = 4 tan af Ge -L)h, Adt+ hy AD] (229 a) 
mg = 4 tan a {hy (A)~ hy (Ame, (xg71) hy (A) 
Haein (OMe Ss - L)E{x -2}h, (A+ by (A)]} (229d) 


maere for 0.5 = AS 1.0, 
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h, A) = 2A - 1 (230@) 
bh, (A) = -A*+ 2A - 3 CET 
HC) SS we ee (230c) 
3 5 6 


This case is identical to that of the unstaggered cascade 
oscillating in phase (d= 0), and the values of pitching 
moment as defined in Eq. (170) and Eq. (171) based on our 
elementary theory for 6 = 0, 8 = 0, and 0.55 AS 1.0 prove 
to be the same as those shown in Eq. (229). 

Drake (1957) extended his work to the case of the airfoil 
oscillating in a wind tunnel with porous walls. Platzer 
(1971) will publish an analysis of this problem using our 
elementary theory. His pitch damping coefficient as defined 


gm Eq. (170), for a slowly oscillating airfoil is, 


C46 = - 4 tan atx? + x (S tan*a-1)+ =(1- tan?) 
eel ae eee mea ce 
*: IO Canc, Wf een) Us Eon at l+otan a? 
ieee @iec an 0) Cot A(1+tan7a))][y]}} (230) 
where 
5 (pak eee 
X (1~A) [ 7% ot ge (At2) 
Y=sx - = (Atl) 
O 2- 


and 0.553 AS1.0. Here 05 6 S © where Go = O represents the 
solid-wall wind tunnel case and ¢ = @ represents a free-jet 
wind tunnel boundary. With these limiting values of porosity 
substituted in Eq. (231) the comparison with Eq. (171b) is 


exact for the case of 6 = 180 and § = 0, respectively. 
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Further the results check with Drake Cig57). “scabilivey boun- 
daries for an airfoil oscillating slowly in a porous-wall 
wind tunnel are given in Figure 7 for a tunnel aspect ratio, 
A = 0.5 and for various values of tunnel porosity. An in- 
crease in tunnel porosity is seen to have a large stabilizing 


influence. 


Me 
15 Ss Unstable 
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Figure 7 


D. AIRFOIL - AIRFOIL INTERFERENCE 

Our elementary solution of the silowly oscillating cascade 
problem can be used to study gther interference problems as 
well. The interference between two airfoils may be of 
interest for certain prgposed space shuttle configurations, 


although two-dimensional flow can only give an indication of 
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trends. Our elementary solution applied to the unstaggered 
cascade oscillating with an interblade phase angle of 180°, 
@a stated before, is equivalent to an airfoil subjected to 
solid-wall wind tunnel interference. If the upper wind tun- 
nel wall is considered a large stationary airfoil (idealized 
as a flat plate at zero angle of attack), the theory can be 
used to determine the non-dimensional pressure distribution 
along the upper surface of a smaller airfoil mounted closely 
below it. Using the non-interference pressure distribution 
over the lower surface of the airfoil, the net pressure dis- 
tribution over the airfoil can be integrated over the chord 
length and the pitch damping coefficient as defined in Eq. 
(171b) can be determined as, 

2 a -~2 tan a a (Ge) = ORS 4GA+ 2 

- + (A?~ 2)(tan?a-1)] + 249-24 - 5(2-A3) 

(tan* a -1)} (232) 


where 


> 
T 


2h cob. eg 


2 
ll 


airfoil distance 
eed 0.5= A= 1.0. 

The stability boundaries for an airfoil undergoing inter- 
ference of this type are shown in Figure 8. The effect of 
interference is strongly destabilizing. This case further 
applies to an airfoil mounted close to only one solid wall 


in a wind tunnel. 
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Figure 8 


E. COMPARISON OF THE METHOD OF CHARACTERISTICS WITH 


ELEMENTARY THEORY 

Values of flow field quantities, 
and imaginary, based on computations 
compared with those quantities using 
as applied to cascades. Comparisons 
values of A, stagger, and interblade 
along both surfaces of the blades as 


field were compared. In the case of 


U, V, and C, both real 
using Program B were 
the elementary solution 
were made for various 
phase angle. Points 
well as in the flow 


Zenomtaterp lade phase 


angle, values of A as low as 0.25 were compared. The restric- 


tion of supersonic leading-edge locus was maintained. All 
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Welues compared were within 2% up to reduced frequencies of 


0.1. Differences between 1% and 2% were obtained only at 


Mach numbers near one (M = ¥1.25) at a reduced frequency of 
O.1. A sample of the flow quantities compared is given in 


fable l. 


bbe C22R G2 C221 C907. 
{} @L.0 uC: SSG Cy M VORCG. S., Sines 





Oo ar Oey, =J 701010 0.060 0.060 
HOO) | —1 29084 =) coo) 0.149 0.150 
3 FOO 1ee5 10 1.500 DOO Ryle l ie 
.00 1536 lees 22 2.189 Del 
00 1.627 1.629 ieatet | 1.162 
0) ee ae ee -2 353 1.633 1.650 
MOO =). 252 es ah I S708 1.738 
#1 if 2 #3 
M= /2 Moe 2 M= /1.25 
A= 0.5 A= 0.5 A= 0.5 
B = 0.0 B = 25.64 B= 13.50 
§ = 36.87 § = 36.87 § = 36.87 
x = 0.0 x = 0.0 x = 0.0 
QO O Deke ie: Ome. es 
Table 1 


KHPrerRoOooog OoCaooaoo © She ee 
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Boriner, comparisons of the out of phase pressure distribu- 
tion on the upper surface of a staggered cascade blade in 
supersonic flow of M = Y1.25 and M = /2 as predicted by the 


method of characteristics and the elementary theory are 


shown in Figure 9 (M= 72, B= .12, aol =moD anid 
Figure 10 (M = 71.25, B = .12, Xo = O, 6 = 0). Here, the 


solid lines indicate the elementary theory, and the dashed 


lines indicate the method of characteristics. 


I fe,/ 8! 
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= Figure 9 


Agreement hetween the two methods is quite good for re- 
duced frequencies up to 0.1 at M = Y2, but this close agree- 
ment is maintained at M = v1.25 only up to a reduced frequency 


of 0.05. Since both methods are based on the same linearized 
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equations of motion, this shows that the influence of fre- 
Guency is greater at the lower supersonic Mach numbers. 

Integrated yalues of lift and moment (Ly and M,) were 
also compared. Here differences as high as 34 were found at 
Mach numbers approaching one; however, for M = /2, the 


values were within 2%. 


ee 
=.2 
5s 
a2 
k=.1 
ae z 
oe k=.05 
=.1 
xa a: i 8, 
Ma k=.0 ee ae 
ere =e a 0 8 x 
=e1l ot 
mig k k=,2 aig 
k=.1 
=,2 
k=,2 


Figure 10 


These results were produced using a grid fineness ratio 
of 100 or one quarter of the capacity of the program in order 
to reduce the amount of computer time required to make the 


computations. 
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F. COMPARISON OF METHOD OF CHARACTERISTICS WITH 
GARRICK-RUBINOW 


Lift and moment calculations obtained from Programs A 
and B at A = 1 were compared with those of Garrick and Rubi- 
now (1946) for both pitch and plunge oscillation. Differ- 
ences were less than 2%. Program B Waa ueed to duplicate 
several different single-degree-of-freedom flutter points 
shown in Garrick and Rubinow's (1946) Figure 22. Again, the 


grid fineness ratio used was 100. 


G. METHOD OF CHARACTERISTICS RESULTS 

The following results were obtained using Program B to 
compute single-degree-of-freedom flutter speeds and frequen- 
cies. Approximately four minutes of CPU time on the IBM-360 
were required to generate one flutter point. A grid fineness 


ratio of 100 was used throughout. 


My, 
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In Figure 11 (Mf = 1.11, A = 0.95, 6 = 0, x) = «5, u = 15.708, 


Tg = ome 0) a typical yarilation of M, with 1/k is shown 
for an unstaggered cascade oscillating in phase with only a 
slight amount of blade interference. Here, the flutter point 
is at a relatively low value of reduced frequency (k = 0.02). 
That this is not the critical flutter condition is shown in 
Figure 12 where the non-dimensional flutter speed is plotted 
versus interblade phase angle for the same parameters as in 
Figure ll. In complete agreement with Lane (1956), the 
critical flutter condition, i.e., the lowest flutter speed, 


for this unstaggered cascade is at an interblade phase angle 


of 180°. 





30 60 90 120 150 180 


Figure 12 
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Figure 13 shows the effect of a small amount of interfer- 
ence on the flutter speed at M = 1.111 when plotted against 


elastic axis position. 


U,, ew 





2 Figure 13 


Again the blade and flow parameters are those as shown in 
Figure 11 and Figure 12. It can be seen that a relatively 
small amount of interference significantly lowers the non- 
dimensional flutter speed over the entire range of elastic 
axis positions tested. 

In Figure 14 (A = .95, 6 = 180, B = 0, np = 500, g = 0, 
rpg = -5) the effect of Mach number on flutter speed is shown 
for various elastic axis positions. An increase in Mach nun- 
ber to M = ¥Y2 significantly lowers the non-dimensional flutter 


speed for a relatively small amount of interference. 
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Figure 14 


In Figure 15 (6 = 180, uw = 15.708, r 


9 5, g = 0) the 
effect of this relatively small amount of interference on 
the flutter speed is again shown, but here at the higher 
Mach number (M = Y2). As before, the non-dimensional flut- 


ter speed is significantly lowered over the entire elastic 


axis range by a small amount of interference. 


Up/ CW, 


a, 





Figure 15 
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Figure 16 (M = /2, A = .95, B= 0, § = 180, x, = 35; 
tg = -9) shows the effect of structural damping on the non- 


dimensional flutter speed for two values of wing density 


parameter. 


UY Cug 
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Figure 16 


With a relatively small amount of interference, the flutter 
speed is greatly increased by a large wing density parameter 
and high structural damping. The high value of wing density 
parameter is more applicable to the cascade case. As defined 
in Garrick and Rubinow (1946), its increase may be interpreted 
as an increase in altitude for a fixed wing density, high 


values indicating supersonic wings at high altitude. For 
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typical supersonic cascades, the wing density parameter is 
of the order of 700. 

The remaining results give variations in non-dimensional 
flutter speed of the staggered cascade (8 = 26.565) at M=/2. 
Results shown are based on a cascade with an interblade phase 
angle of 180°. It must be re-emphasized that this interblade 
phase angle is not necessarily the one that gives the criti- 
cal flutter condition; however, several points at random were 
chosen and the critical flutter condition found for each. 

All were within 1% of flutter speed at 6 = 180 and within 
10° of 6 = 180. This was not done for all points shown in 


an effort to reduce the total amount of computer time used. 
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Figure 17 


Dnetieure 17 (i = 500, ¢ = GZ. = .5) the effect of an 


*9 
increasing amount of interference on flutter speed is shown 


for various elastic axis positions. The cascade is assumed 
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to have a high wing density parameter and 24 structural damp- 
ing coefficient. An increase in the amount of interference 


is seen to lower the non-dimensional flutter speed. 
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Figure 18 


Figure 18 (A = .6, es -5) shows the effect on the non- 
dimensional flutter speed of increasing the wing density 
parameter then the structural damping for various elastic 
axis positions. Without structural damping, an increase in 
wing density parameter had little effect on the flutter 
speed. With 2% structural damping, however, the flutter 
speed was greatly increased. The effect of structural damp- 
ing is shown further in Figure 19 (A = .6, Xo * a" T= ioe 
Mere, the variation of Elutter speed with structural damping 
is shown for three yalues of wing density parameter. As 


can be seen, increasing the structural damping increases the 
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Figure 19 


difference in flutter speeds for the three values of wing 
density parameter. Thus, the accurate determination of the 
wing density parameter becomes increasingly more important 


with increased structural damping. 
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Figure 20 
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Figures 20 and 21 show the same Parameters as Figures 18 
and 19, respectively, but here, there is more interference 
(A = .25). The results are relatively the same as those 
Shown in Figures 18 and 19. Flutter speeds overall are lower 
due to the increase in the amount of interference, and the 


effect of structural damping is more pronounced. 
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“VLIL. CONCLUSIONS AND RECOMMENDATIONS 


The elementary theory for cascades with supersonic lead- 
ing-edge locus is valid for low frequencies and can be used 
to give rapid calculations of lift forces and pitching mom- 
ents of cascades as well as their torsional stability boun- 
daries. These results can immediately be applied to a wide 
variety of wind tunnel and airfoil interference problems. 

The method of characteristics can be used in conjunction 
with the high-speed computer to give rapid, accurate flutter 
calculations of cascades with supersonic leading-edge locus 
oscillating at arbitrary frequency. This applies to torsion 
and bending flutter calculations as well as to torsion flut- 
ter alone, for the only diffexrences in the calculations are 
those due to different values of normal velocity at the blade 
surface. In each case, torsion or bending, the method of 
calculating the flow field and determining the lift forces 
and moments is identical. 

Results of single-degree-of-freedom (torsional) flutter 
calculations show that interference on cascade blades due to 
reflected Mach waves from adjacent blades lowers flutter 
speeds significantly and in proportion to the amount of this 
interference. The presence of structural damping, however, 
can partially alleviate this adyerse situation by increasing 


Flutter speeds, provided the wing density parameter is large. 
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Recommendations for future study include the development 
of the elementary theory for bending oscillations, and the 
generation of numerical results for two-degree-of-freedom 
flutter using Program A. O£f primary importance, however, is 
the development of a method of characteristics approach for 
cascades with subsonic leading-edge locus. This is the criti- 
cal problem today in the design of the supersonic compressor. 
With this problem properly formulated, the method of charac- 
teristics should be able to give rapid, accurate flutter 
predictions for cascades oscillating in this complicated 


flow field. 
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APPENDIX A 


COMPUTER PROGRAM VARIABLES 


L. LOGIC VARIABLES 

The following variables and their definitions, in alpha- 
betical order appear in the logic statements of the programs. 
The asterisk is placed behind those variables that appear 


only in Program A, 


FCONST* Denotes type of oscillation: TRUE - plunge 
mode, FALSE - pitch mode. 

LATNWL When computation reaches upper airfoil, it 
takes on the value "1", otherwise, "0". 

LCO Takes on the value "1" when lower airfoil 
computation is complete, otherwise "0". 

ICOUNT Incremented only when MACHLN is called by 
being set equal to JLINE, zeroed when HIFOIL 
is called. 

IFIN Seteequal to, 1) whemeinitecraGion, on lower 
airfoil is complete, otherwise "0". 

LHAVEP The number of the computation in any row 
Star tinesiesen “Or: 

LJUNC The grid point at the intersection of the 
initial right and left-running Mach lines. 

LJUNT Equal to IJUNC, needed to preserve the value 
initially calculated in case of multiple 
flow field calculations. 

ILINE The number of grid points in the row (right- 
running Mach line) being computed, plus one. 

IMAGRT* If X is negative, set equal to "1", other- 
wise "O", 

LREF Set equal to "1" if lower airfoil point is 
in Zone I, otherwise "0". 

Per bit Set equal to 1/2 of the number of grid points 


minus the number of stagger points. 
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ISWICH 


LSWICH 


ITEM 


IZOT* 


JCOUNT 


JFIN 


JLINE 


JSWICH 


JSWICH 


JTEM 


KCOUNT 


KOUNT 


LCOUNT 


MCOUNT 


Semmcntmais to b' soe "2" alternately starting 
with "1" and changing when upper airfoil cal- 
CUbattomecoes intomasdttnerent zone. The 
upper airfoil zone marker. 


Switch yariable used to choose the row being 
calculated, alternately set equal to "1" or 
"2". Used as a second subscript. 


Loop variable for pitch and plunge mode cal- 
Cuda t Lon, 


Scie CUaent Omen © oC nteoweErequency rats o 
is imaginary, otherwise "0". 


Set equal to "0" when in Zone I, incremented 
when MACHLN is called in an interference 
region. 


Set equal to "1" when integration along 
upper airfoil is complete, otherwise "0". 


The number of grid points in the row being 
computed. 


The same as ISWICH but used in LOFOIL. The 
lower airfoil Zone marker. 


Switch variable used to choose the preceding 
line that was calculated, set equal to the 
Opposite value as that of ISWTCH. Used as 

a second subscript. 


Loop variable for frequency range increment- 
ing. 


Used in the interference calculation in 
MACHLN, set equal to "JCOUNT-1". 


Set equal to "0" until HIFOIL is called, then 
set equal to the number of grid points minus 


one. 


The number starting with "0" of each point 
in a zone of interference in LOFOILL. 


Same as LCOUNT but used in HIFOIL. 
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IL. QUANTITY VARIABLES 


The following variables in the programs take on the 


values as defined below. All have real values unless other- 


wise stated. The asterisk is placed behind those variables 


that appear only in Program A and the accent is placed behind 


those variables that appear only in Program B. The dimen- 


Eronearvariables are listed tirst. 


W22R. U22T, 


Usenmeu331, 


ee YeCa Oi 2 ) 


AI 
BI 
CAPX' 


CPi 


CPRI 


CPI2 


CPR2 
DAMP ' 


DELCI 


V22R, V221, C22RF7° E2295 (400 53). The velo- 
cities, real and imaginary, of a point in 
the program. They are equivalent to the 
respective Teipel amplitude functions, 
Usor through Coot’ 

V3orenV 331,51C35r. C331(200,2). Same as 
above, but used to calculate the flow field 
in Zone I' when the phase angle is non-zero. 


The x and y coordinates of a point based on 
the chord length. 


A. as defined in Eq. (96a). 


if 
By as defined in Eq. (96b). 
X as defined in the Table of Symbols. 
The imaginary part of 7 (eu amt re the 
preceding point calculated on the 
upper airfoil. 
-id : 
The real part of e (Co 5ptil, oy) of the 


preceding point calcu- 
lated on the upper airfoil. 


The imaginary part of et iG yottc ae the 
point being calculated on the 
upper ata om ss. 


The real part of e *? eopgatee iC,,,Jof the point 
being calculated on the upper airfoil. 


The total system damping as defined in Eq. 
C2 ear 


The value of Coot at x = l. 
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DELCPI 


DELCPR 


DELCR 


DELL1, DELL2, DELM1, DELM2. 


DELTA 


DELTAS 


DISC# 
Dots ER 
FACTOR' 


FACTOR* 


FAZE 
FIN 
FMANGL 


FNGDPT 
FNGRDN 


FSTRMN 


ORs OS 
GSUBA* 
HAFREQ* 


HDSTRL 


The imaginary part of e+? @ELCR TIDEEEL) on 
the upper airfoil. 


14 


The real part of e ~*(DELCR + iDELCI) on the 


upper airfoil. 


The value of C ate =). 


22R 

The change in lift (real 
and imaginary) and moment, respectively, due 
to extrapolation to x = l. 


The phase angle in radians. 


The size of the length increment along a 
Mach line. 


The discriminant of Eq. (206). 
The size of the Ax increment. 
Set equal to Ax/k?. 
In the MAIN, set equal to Ax/k? for the 
pitch mode, and set equal to 1/2 Ax/k? for 
the plunge mode. In FLUTTER, set equal to 
DilsCe= 1C>:. 

R 
The phase angle in degrees, 
The chord length. 


The Mach angle in radians. 


The real number value of the number of grid 
points on the initial Mach line. 


The real number value of the grid fineness 
ratio, equal to FNGDPT. 


The freestream Mach number. 


Intermediate factors used in velocity calcu- 
lations in MACHLN, LOFOIL, GENFPT, and HIFOIL. 


The structural damping coefficient, g, as 
defined in the Table of Symbols. 


The ratio of the natupaterrequency Of sthe 
alreoll in plunge esto mth atm piste. 


The increment, Ax/2. 
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INCRE The amount that 2/k is incremented each time 
the lift and moment is calculated. 

KSPR... »KOGL Kyope tos 6y as defined in Eq. (95). 

iy) 2 The real and imaginary lift coefficients, 
respectively. They become Ly and L, or L, 
and Ly as defined in Eq. 
(186)and Eq. (187) depending on the type of 
oscillation. 

L3, ba* The real and imaginary lift coefficients, 
respectively, due to pitching oscillations. 

M1, M2 The real and imaginary moment coefficients, 
respectively. They become My and M9 or M. 
and M, as defined in Eq. 
(188)and Eq. (189) depending on the type of 
oscillation. 

M3, M4% The real and imaginary moment coefficients, 
respectively, due to pitching oscillations. 

MUU The wing density parameter as defined in the 
Tabte of symbpows., 

NGDPTS The grid fineness ratio plus one. 

NGRDFN The grid fineness ratio, it must be an even 
number less than 400. It is an integer. 

NSTPTS The number of stagger grid points, an integer, 
set equal to B/Ax plus one. 

OMEGA The non-dimensional flutter frequency as 
defined in Eq. (208)- 

OMEGAA 2 as defined in the Table of Symbols. 

OMEGAH* 7 as defined in the Table of Symbols. 

REDFRQ ihre reduced frequency, 

RFREQ PEs Le 

ROOTX The square root of X. 

ROOTX1* The square root of themedrst root of Eq. (206). 

ROOTXK2* The square root of the second root of Eq. 
CZUG0)r. , 

ROOTXI* The square root of 9a. (207). 
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RSUBA The radius of gyration, tg» as defined in 
the Table of Symbols. 


RTOSPH Thestatio of specitic heats. 
S The square root of M2-1, 
STGANG The stagger angle in degrees. 
STGR The amount of stagger, B. 

T M2~1, 


TEST1, TEST2* The difference between the respective 
FOOtS Of EQ ec UO) and chat of whe. “(20/7 


TEST3, YEST4* The values of TEST1 and TEST2, respective- 
ly, for the previous frequency computation. 


TNWDST The perpendicular distance between the two 
air ois sed. 





ZOECRD The chord length of the top airfoil, x-B. 
TRNGLH The increment, =Ax Gos! 40. 
M2 
U The factor, cos(k x). 
M?7-1 
UF The non-dimensional flutter speed as defined 


me de, MC LUG. 


Um UL, VR, VI, CR, Cl. Mhe aatues of Vv 
pectively, for calculations 
on the initial Mach line between the point of 
intersection of the initial left and right- 
running Mach lines and the upper airfoil. 


lin??? *? yr res- 





2 
V The factor, “sim eek zZ xa 
M?-1 
VIPANL The imaginary part of the normal flow velo- 
city at thé advioilecurtace. 
VRPANL The real part of the normal flow velocity at 
tie ate foes Ure caer, 
W Mor 
XI : Xy as defined in Eq. (207). 
XLENGTH The length of the initial Mach line. 


XR1, XR2* The roots of Eq. (206). 


go 





XSUBA*# The distance from the elastic axis of the 
alrtoll tometsmcenter of gravity. 


XSUBO The distance from the leading-edge of the 
airfoil to its elastic axis. 
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APPENDIX B 


FLOW DIAGRAM - PROGRAM A 


I. MAIN PROGRAM 


Read/Write DATE 
CALL INPUT 





IZOT = 0 


100 
KTEM = 1,4 
1005 
JTEM = 1,50 


F CONST FALSE. 


1002 
ITEM = 1,2 


IcO = 0 
KOUNT = 0 
CALL INITIAL 


1000 
CALL COMPXY 





TATE 0.0 


o AND. 
LATNWL. EQ. 






CALL MACHLN 


THAVEP. EQ.0 
. AND. 
LTATNWL. NE. O 






CALL HIFOIL 


TATNWL. EQ. 1 


~ AND. 
JHAVEP, EQ. KOUNT 


LTHAVEP. EQ. JLINE 


ig T 


. | CALL LOFOIL 






oF 


















CALL GENFPT Leal 


IATNWL. EQ. 1 
. AND. 
HAVEP, EQ. KOUNT 
F 


LTHAVEP.NE. JLINE 
F “i 


(X (LHAVEP+1, ISWICH)+DSTSTR) .GT.FIN 
ILINE.EQ.NGDPTS 


JLINE-1.LT. ITJUNC 















TATNWL = 1 
JLINE = NGDPTS 
KOUNT = NGDPTS-1 


ILINE.EQ.NGDPTS Mis. oc = 1 


ILINE = ILINE+1 
JLINE = JLINE+1 








ISWICH. EO. 1 








THAVEP = ITHAVEP+t1 


ICO.EQ.1 
-NOT.FCONST TNC. £0. 0 
F 
Complete Complete ICO. EQ.1 
eee enh L3,L4,M3,M4 _AND 
T IHAVEP. GT. LJUNC 
F 
FCONST = .TRUE. 


1002 
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CALL FLUTER 








JTEM.NE. 1 
IMAGRT. NE. 1 OR. 
IZOT.NE. 1 


IZOT = 
Define TEST3, 
TEST4 












RFREQ = RFREQFINCRE 
REDFRQ = 2/RFREQ 


ABS (TEST1).LE. LE-05 
Write OMEGA, UF, JTEM 7 OR, 
XSUBO = XSUBO + 0.1 ABS (TEST2) . LE. 1E-05 
F 


Define TESTI, 
feor 2 










1100 TEST1.LT.-1E-05 
F 

STOP TEST3.LT.-1E-05 

END 


TEST2.LT.-LE-05 











TEST3.GT.1LE=-05 


TEST4.LT.-1LE-05 
F T 


TEST4.GT.1E-05 


RFREQ = 
INCRE = 


= INCRE.LT.1E-05 


RFREQ-1 
INCRE/5 
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II. SUBROUTINE INPUT 


Read/Write NAMELIST 
Calculate DSTSTR, NSTPIS 


STGR-DTSTSR.GE.0.50 | NSTPTS = NSTPTS+1 
NGDPTS /2.EQ.NSTPTS-1 NSTPTS = NSTPTS-1 


NGDPTS /2.EQ.NSTPTS+1 NSTPTS = NSTPTS+1 





Calculate STGR, STGANG, IJUNT,NSTPTS, 
ISPLIT, S, DELTA, OMEGAA, OMEGAH 
Write NAMELIST 

RETURN 

END 








III. SUBROUTINE INTIAL 


Define IJUNC, AI, BI 
~ NOT. FCONST 


Define VRPANL, VIPANL 


Calewlate U22R> C227 (1,1). XC tay ie) 







Define VRPANL, VIPANL 


ISPLIT.NE.-1 


DELTA. LE. LE-O5 







Calculate L1,L2, M1,M2 


Callenwate LE.L2. MieMZ 


100 


a — ie 
- ; 
— 
a 
: 





DELTA. LE, LE-05 


Je 


. NOT. FCONST 


ic 







Define VRPANL, VIPANL 


Calculate U33R - C331(1,1) 
ISPLIT.EQ.-1 


Define VRPANL, VIPANL 







GCarecurace Lt,L2, MI.M2 


Hof 


rm rb 


ore iii 





101 





ITV. SUBROUTINE MACHLN 


ICOUNT = JLINE 
7 ICOUNT. LT. TJUNC 


F 
JCOUNT = JCOUNT+1 
KCOUNT = KCOUNT-1 


KCOUNT. NE. 0 







Define UR,UI, VR, VI, CR, CI Define UR,UI, VR, VI, CR, CI 
Calculate U22R ~ C22I(1,ISWICH) 
Define W,T 


» NOT. FCONST 





Define VRPANL, VIPANL 
Calcumate UzZR > C221(1. ISWICH 


Define VRPANL, VIPANL 


COUNT. NE. ISPLIT DELTA. LE. LE-05 
OR. OR. 
DELTA. GT. LE-05 T JLINE. GT. ISPLIT T 
F F 
Define Y22R > C22I(1.3 = » NOT. FCONST 
F 
Define VRPANL, VIPANL Define VRPANL, VIPANL 


ICOUNT NE. ISPLIT Calculate U33R ~ C331I(1, ISWICH) 





Define U22R —- C221I(1,3 
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V. SUBROUTINE HIFOIL 


Define TOPCRD 
JLINE. EQ.NGDPTS-1 ' wa LATNWL = 0 


ICOUNT. EQ. 0 





JCOUNT.NE.O 





Define V22R,V221I (1, ISWICH) 


Define K12R,K121I,J 






Calculate U22R, U22I, C22R, C221(1,ISWICH), 
CPR1, CPI1, CPR2, CPI2 


Jin EQ. 1 


COUNT. EQ.0 
«AND. 


OPCRDEDSTSTR. GT.1.0 T 
F 


LCOUNT. EQ. 0 
OR. 


OFCEDTDSIS—IR. GI.1.0 








Calculate L1,L2,M1,M2 


Cateulate i1.L2.Ml] M2 


LCOUNT = LCOUNT+1 
130 
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130 


ISWICH.NE.1 














LCOUNT. NE. JCOUNT 


LCOUNT. EQ. LJUNT 


Calculate L1, L2, M1,M2 
LCOUNT = 0 


ISWICH. EQ. 






ISWICH = 2 


ISWICH = 


OECRDzDoISIR. LE. 1.0 


Catculote. Ul, LZ, Ml, M2 
JFIN = 1 


ICOUNT = 0 
RETURN 
END 
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VI. SUBROUTINE GENFPT 


IATNWL. EQ.0 







Define K12R, K12I, K34R, K34I Define K12R, K12I, K34R, K34I1 


Calculate U22R — C22I(1I, ISWICH) 


COUNT. NE. ISPLIT 
ORE 
DELTA. GT. LE-05 












Refine U2Z2R = C221(1,3) 


DELTA, LE. LE-05 
~OR. 


JLIOVE. GT. [SPLIT 





Define K12R, K12I, K34R, K34I, K56R, K561 
Calculate U33R — C331(1, ISWICH) 


JLINE.NE.ISPLIT 
DetimnesU22R —.. C221 (1,3) 


UO), 





VII. SUBROUTINE LOFOIL 


- NOT. FCONST 
Detime VZ2R,V221 (LINE, ESWweEen) 


Define K12R, K12I, K34R, K34I 
LATNWL. EQ.0 







Define V22R,V22I(ILINE, ISWICH) 






Define K56R, K56I Define K56R, K56I 


Calculate U22R,U221, C22R, C221 (ILINE, ISWICH 


LCOUNT.NE. ISPLIT 


DELTA. GT. 1LE-05 







Define U22R ~ C22T(ILINE, 3 





DELTA. LE. LE-05 
OR, 
i JLINE. GI. ISPLIT 


NOT. FCONST 







Define V33R, V33I (ILINE, ISWICH) 


Cabeulate U33R,U331,C33R,C331n( LLINESESWICH) 
JOINE SNES ESPLIT 


Derane U2Z2R = C220 ( TLINE, 3) 


[88 


Define V33R,V33I(ILINE, LSWICH) 







106 


: 7 





Define CPR1, CPI1, CPR2,CPI2 


IFIN. EQ.1 
aly 


F 


Calculate L1,L2, Ml, WZ 


JLINE. NE. ISPLIT 
. AN 








X(K, L)+DSTSTR.LE.1.0 


Catculate LI,LZ2, ML, MZ 
IFIN. EQ. 1 


MCOUNT. EQ. 0 
. AND. 

X (CILINE, ISWICH)+DSTSTR.GT.1.0 ak 

F 

140 

MCOUNT. EQ.0 

JOR. 
X(ILINE, ISWECH)+DSTSTR.GT.1.0 

Galeulate LIZuZ ML. M2 

EREES EO! 


X(LGINE, ISWieGH EDSTsirn. Gr.l.c 








GCalculaceri. L2. MiFiZ 













MCOUNT = MCOUNT+1 


JSWICH, NE, 1 







MCOUNT.NE. ISPLIT+1 


TNE, : 
F 2h 
GREF, EQ) 
T F 
120 


Bu(0) 7 





MCOUNT = 0 
JSWICH. EQ. 1 


JSWICH = ISPLIT. EQ. -1 7 
|JSWICH = 1} 


ISPLIT. EQ.~1 130 JSWICH 





JSWICH = 2 














JCOUNT. NE. 0 
aOR. 
DELTA. GTI.1E-05 





NeEPrs. EQ. 1 


F 
X(ILINE, ISWICH)+DSTSTR.GT.1.0 = 195) 
F 
MCOUNT. EQ. 0 —>{110] 
T 
X (ILINE, ISWICH)+DSTSTR. CT. 1.0 


MOURNE toWECH Do tSih., LE ano . 140 
F 
Calculate DELL1, DELL2,DELM1,DELM2 


Caleulace Ll k2. ML, M2 


Calculate L1,L2, Ml, M2 
IFIN = 1 


URN 


DELTA. GT.1E-05 
TREF.EQ.0 


MCOUNT. GT. ISPLIT+1 RET 
END 
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VIII. SUBROUTINE COMPXY 


THAVEP. EQ. 0 















Define I 
alculate X,Y(1I, ISWICH) 


Define I 
Calculate X,Y(1I,ISWICH 


LATNWL. EQ. 1 
AND. 
JLINE.NE.NGDPTS-1 
19 


Define I 
Calculate X,Y(1I, ISWICH) 










IX. SUBROUTINE FLUTER 


IMAGRT = 0 
Calculate CR,CI 


OMEGAH. LT, 1E-06 
20 


R. T 
OMEGAA. LT. LE-06 
F 
Calculate DISC Calculate XR1,XR2 


DISC, LE. -1LE-06 










Calculate XR1,XR2 
Calenlate. x1 


XR1.LE.-1LE-06 
~OR. 


XR2.LE.-1LE-06 
~OR, 
XI. LE.-1E-06 


Calculate OMEGA, UF 
Write RFREQ, ROOTX1, ROOTX2,ROOTXI 





IMAGRT = 1 
Write RFREQ 
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